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CALCULATION BY A FINITE-DIFFERENCE METHOD OF 


SUPERSONIC TURBULENT BOUNDARY LAYERS 
WITH TANGENTIAL SLOT INJECTION 

By Ivan E. Beckwith and Dennis M. Bushnell 
Langley Research Center 

SUMMARY 

A method has been developed for calculating compressible turbulent boundary 
layers with tangential slot injection of homogeneous gas species. The partial differential 
equations for the mean motion are solved by an implicit finite-difference method. The 
turbulent-flux terms are modeled by means of eddy-diffusivity and mixing-length con- 
cepts. The magnitude and distribution of the mixing length across the boundary layer 
are determined from the computed characteristics of the free-mixing region between the 
injected jet and the initial boundary layer. The development of this mixing region and its 
interaction with the wall and external boundary layer are calculated from the species- 
conservation equation which, in the present computer program, is used to calculate the 
behavior of a trace species. 

The method is described in detail, and the appropriate initial and boundary condi- 
tions used to calculate the concentration of the trace species and the corresponding 
mixing-length relations are also described. The numerical-solution procedure is simi- 
lar to that contained in a previous computer program developed at NASA Langley Research 
Center to solve conventional boundary -layer problems. The modifications to the previous 
program required to solve the present slot-injection problem are described and listed in 
FORTRAN IV in the appendix of this report. 

Comparisons of predicted velocity profiles, boundary-layer thicknesses, heat trans- 
fer, skin friction, and recovery temperatures (or effectiveness) have been made with 
results from four previous experimental investigations at free-stream Mach numbers of 
3 and 6. Agreement with experimental data was generally good. In particular, good pre- 
dictions were obtained of velocity profiles and surface properties in the region just down- 
stream of the slot where no previous method has been successful even for low-speed 
flows. The generally good agreement obtained between theoretical results and experimen- 
tal data indicates that reliable predictions can be obtained for slot-injection flows by the 
present method even when the inherent restriction of small normal pressure gradients in 
boundary -layer theory is violated. 



INTRODUCTION 


Injection from slots into turbulent boundary layers has long been advocated as a 
method of controlling the downstream -wall temperature and skin friction. The injection 
of heated air is a well-known and effective deicing method. (See ref. 1.) In hot environ- 
ments, the injection of cool air from either a slot or porous strip is usually referred to 
as film cooling since it provides thermal protection for the downstream surface. (See 
refs. 2 and 3, for example.) When the injected fluid is directed downstream through a 
rearward-facing tangential slot, the skin friction generally increases or decreases 
according to whether the specific momentum of the jet is greater than or less than that 
of the free stream. With jet momentum greater than free-stream momentum, the flow 
configuration is often referred to as a wall jet, and the resulting increase in skin friction 
may be utilized to delay separation. A recent review of experimental data and integral 
theories for such flows is given in reference 4. The cooling effects of wall jets have 
been investigated for incompressible flow in reference 5 where it was shown that 
Colburn's Reynolds analogy, expressed in terms of the maximum velocity in the wall jet 
and applied to empirical correlations of skin-friction data, gave satisfactory predictions 
of the measured heating for this type of flow in the region far downstream of the slot. 

The flows to be considered in the present report will be primarily those where the 
momentum of the tangentially injected fluid is less than that of the free stream, and some 
reduction in skin friction might be expected. Because of the possibility of simultaneous 
thermal protection and skin-friction reduction, there is continuing interest in this type of 
film cooling and, consequently, a large amount of literature is available on the subject. 

Typical experimental results and correlations of film cooling effectiveness for 
incompressible flows with constant free-stream velocities are given in references 1 and 
6 to 9. Comparisons in reference 8 of data from some of these investigations and the 
subsonic compressible data of reference 2 indicated the effectiveness (a normalized form 
of the local adiabatic-wall temperature) may vary considerably with changes in slot con- 
figuration, ratios of slot height to initial boundary-layer thickness, mass -injection ratios, 
and compressibility effects. These comparisons have led to several investigations 
(refs. 10 to 13) concerned with various aspects of slot geometry, such as slot width or 
height, injection angle, and lip thickness. In practical applications of film cooling, sev- 
eral other effects must also be considered, such as those caused by foreign-gas injection, 
variable free-stream velocities, and multiple slots. 

The influence of foreign-gas injection, which may introduce large density variations 
even in low-speed flow, has been investigated in references 13 and 14. The effects of 
variable free-stream velocity have been reported in references 15 and 16, where it is 
indicated that film cooling effectiveness is reduced by strong accelerations downstream 
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of the slot. Experimental data and approximate methods for computing wall temperature 
downstream of multiple slots are presented in references 17 to 20. 

Most of the references for slot injection discussed thus far are for low-speed flows, 
where compressibility effects due to high velocities are negligible. Experimental data 
for higher speeds but subsonic free-stream Mach numbers are reported in references 2 
and 16. In reference 20, injection into supersonic turbulent boundary layers has been 
treated by reference-enthalpy methods applied to correlations of low-speed data. Exper- 
imental data and correlations or approximate theoretical analyses for tangential slot 
injection into supersonic turbulent flows are also available in references 21 to 29. 

Detailed boundary-layer-profile data for these conditions are available only in refer- 
ences 21 to 23, 

Velocity-profile data for injection from a sonic slot into a flat-plate turbulent 
boundary layer at a free-stream Mach number of 3.0 are given in reference 21. The 
main purpose of this investigation was to measure skin friction and drag, and no data 
were obtained on film cooling effectiveness since the total temperatures of the slot and 
free-stream flow were nearly equal (294° K (530° R) and 317° K (570° R), respectively). 

The investigations of references 22 to 24 were at a free-stream Mach number of 6.0 
with free-stream and slot-flow total temperatures of about 444° K (800° R) and 294° K 
(530° R), respectively. However, adiabatic temperatures were not measured directly in 
any of these investigations; instead, the film cooling effectiveness was inferred from 
measured heating rates and calculated heat-transfer coefficients from the flat-plate 
reference-enthalpy method of Eckert (ref. 30). Comparisons in reference 22 of these 
results at Mach 6 with previous data indicated that film cooling may be much more effec- 
tive in hypersonic flow than in subsonic flow. However, the indirect procedure for cal- 
culating adiabatic-wall temperatures, as used in references 22 to 24, is questionable 
because the calculated heat-transfer coefficients based on flat-plate equilibrium flow do 
not apply to slot-injection flows. The actual values of heat-transfer coefficients in a 
boundary layer perturbed by rapid changes in surface temperature (such as occurs with 
slot injection) are known to deviate considerably from the values for an undisturbed 
boundary layer as shown, for example, by the results of reference 31 and by the analysis 
of reference 32. The analysis of reference 32 was used to predict the adiabatic-wall- 
temperature distribution for the experimental investigation of reference 1 by means of a 
specified step function in surface heat transfer. This step function simulated the heat 
addition at the slot. While the predicted values of effectiveness were about 20 percent too 
high, the predicted trend with downstream distance and mass -injection ratio was the same 
as that of the data in the far-slot region. Hence, the dominant factor that determines the 
adiabatic -wall temperature in the far-slot region for low-speed flows is the initial energy 
content of the flow near the wall. However, in the near-slot region the momentum transfer 
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between the jet and initial boundary layer is probably of crucial importance. In particu- 
lar, Colburn's Reynolds analogy factor may not be applicable in this region because of 
the different structure and rapid changes in the velocity and thermal mixing layers. 

Also, when the injected momentum is less than the free-stream value, in contrast with 
the wall jet (ref. 5), there is no consistent reference velocity close to the wall that can be 
used to formulate a simple Reynolds analogy expression. 

Adiabatic -wall temperatures downstream of an injection slot were measured 
directly for the first time in hypersonic flow in the investigation of reference 29 which 
was also at a free-stream Mach number of 6.0. These results indicated that significant 
improvements over subsonic and low Mach number supersonic results can be expected in 
both the extent of the thermal protection and in the relaxation rate to undisturbed equilib- 
rium temperatures. 

It is not known whether these increased efficiencies of slot cooling observed in ref- 
erences 22 and 29 are due to slower mixing rates in the hypersonic boundary layers or to 
the relatively small ratios of slot height to initial boundary-layer thickness (0.01 to 0.2) 
used in these investigations. One approach to such problems is by means of finite- 
difference solutions of the flow field with the complete initial and boundary conditions 
included. At the present time, solution techniques for the complete equations of motion, 
including the normal momentum equation (see ref. 33, for example), are impractical for 
engineering applications because of the lengthy relaxation procedure and large number of 
mesh points required to obtain good accuracy within an entire flow field for which all 
boundary conditions must be specified. Furthermore, the results would be only as good 
as the models of the turbulent-flux terms. 

When the normal pressure gradients can be neglected, the general equations reduce 
to the conventional boundary-layer equations, which are parabolic and hence can be solved 
by forward-marching techniques with only the initial conditions and boundary conditions in 
the free stream and at the wall required. Models for the turbulent-flux terms can then be 
developed and tested with these simpler flows and computing procedures. Both explicit 
(ref. 34) and implicit methods (refs. 35 to 37) have been applied successfully to turbulent- 
boundary-layer flows without slot injection. 

The method of reference 37 has also been used extensively to calculate low-speed 
turbulent-boundary-layer flows with tangential slot injection. (See refs. 14 and 38, for 
example.) A simple two-step ramp function for the Prandtl mixing length was used in 
these calculations, and good predictions were obtained for profiles and film cooling effec- 
tiveness except in regions less than 20 slot heights downstream of the slot. Experimental 
data for low-speed flows (ref. 39) show that normal pressure gradients are small for 
downstream distances from the slot of greater than 20 slot heights when the ratios of 
injection velocity to free-stream velocity are less than unity. 
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For slot injection into supersonic flow, however, large disturbances which propagate 
far downstream are generated at the slot (refs. 25 and 27) unless the slot static pressure 
is carefully matched to the local-stream static pressure, as in reference 26. As pointed 
out in reference 27, the flow field that results when large differences in slot and stream 
static pressure are present cannot be computed with the conventional boundary-layer 
equations. The shock and expansion waves that dominate the flow for this "mismatched" 
pressure condition can only be accurately accounted for by the complete equations of 
motion. 

The purpose of the present investigation is to determine if the characteristics of 
supersonic turbulent boundary layers with tangential slot injection can be predicted with 
the finite-difference method of reference 36. Since the normal-momentum equation is 
neglected, the results may not apply when the static pressure at the slot exit differs from 
the local-stream static pressure, or when the slot -injection flow is inclined at an appre- 
ciable angle to the local-flow direction. In an attempt to improve predictions in the 
important region just downstream of the slot, the simple-ramp mixing-length distribution, 
as used in references 14, 37, and 38, has been modified to account for the rapidly changing 
flow structure in this region. Again, the results may be questionable when the slot-lip 
thickness is appreciable compared with the slot height since for incompressible flow the 
effects of the thick lip influence the flow far downstream to distances of 100 or more slot 
heights. (See refs. 12 and 13.) The present theoretical predictions will be compared 
with the experimental data of references 21 to 24 and 29 for free-stream Mach numbers 
of 3.0 and 6.0. While the present method has not yet been applied extensively to subsonic 
or incompressible slot-injection flows, preliminary results for such flows (not reported 
herein) indicate that reliable predictions can be obtained. For the present applications, 
it should be noted that in the investigations of references 21 to 24 and 29, the restrictions 
mentioned previously regarding matched static pressures, small injection angles, and lip 
thickness were generally violated to some extent. 

Modifications to the original boundary -layer computer program (ref. 41) reqioired 
to compute the present slot-injection flows are given in the appendix by Barbara A. Hixon 
and Dennis M. Bushnell. 


SYMBOLS 

Measurements and calculations were made in the U.S. Customary Units. They are 
presented herein in the International System of Units (SI) with the equivalent values given 
parenthetically in the U.S. Customary Units. 

AdP 

A Van Driest* s damping parameter, — (see eq. (3)) 

pVWp 
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damping constant, 26 for zero mass transfer by normal injection at wall 
ratio of mixing length to width of various flow regions (see eqs. (14) and (18)) 

2Tw 

skin-friction coefficient, 

o 

*^e e 

^i 

mass concentration of species i, — 

specific heat at constant pressure 
diffusion coefficient of species i 

velocity -profile variable, — 

^e 


mixing-length function from reference 36 (see eqs. (4)) 

C, - C. 


normalized concentration profile, 


'1 ^i,w 
^i,e “ ^i,w 


'^i 

ratio of local to edge concentration of species i, — 

^i,e 


total enthalpy, h + — 


static enthalpy 

constant in Prandtl's mixing-length relation, taken as 0.4 herein (see eq. (4a)) 
reference length 
mixing length 


Mach number 



m 


molecular weight 


NLe 

Npr 

%e 

Nsc 

Nst 

P 

q 

R 

r 

s 

T 

t 

u,v 


Lewis number 
Prandtl number 

PgUgX 

Reynolds number, 

Me 

Schmidt number 
Stanton number 
pressure 

heat-transfer rate 

universal gas constant 

radius from axis of symmetry 

slot height (see fig. 1) 

absolute temperature 

thickness of slot lip (see fig, 1) 

velocity components in x,y directions 


V = V 


+ 




p 


W height of concentration mixing region 

Ay m Q 

X correlating parameter for streamwise distance and mass flow, — 

s 

x,y boundary -layer coordinates, parallel (in flow direction) and normal to surface. 

For free turbulent flows, y is distance from plane of symmetry for jets 
and distance from plane where u = + ^min) for half jets. 
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ratio of specific heats 


boundary-layer thickness evaluated where F = 0.995 unless otherwise noted 
initial total flow thickness, 6 q + t + s 


6 

(1 - F)dy 


For the present problem, the lower limit is replaced by y at 
G = 0.01 to give 6in,s) 


eddy viscosity (see eq. (1)) 


ratio of local stagnation temperature to free-stream stagnation temperature, 

Tt/Tl,e 


recovery-temperature effectiveness. 


'^t,e ~ Taw 

T - ^ 

^t,e ^t,j,o 


total enthalpy variable. 


He-Hw 


1^(1 - F)dy 


F(1 - F)dy (Por the present problem, the lower limit is replaced by y at 
G = 0.01 to give %n,s) 

mixing-length factor for pipe or channel flows 

specific mass-flow ratio of injected flow to free-stream flow. 



molecular viscosity 


density 


T 


shear stress 


(jj recovery factor (see eq. (32)) 

Subscripts: 

A air 

av average or mean 

aw adiabatic wall 

b outer part of boundary-layer region with slot injection 

c concentration 

e "edge” of boundary layer where specified boundary-layer profile parameters 

F, 0, or g -* 0.9999 

f far-wall region of boundary layer 

i species 

j jet region 

m free turbulent mixing region 

max maximum 

min minimum 

n near -wall region of boundary layer 

o initial values or station, usually at slot exit 

r reference 

s slot injection flow 


9 



T turbulent 

t local stagnation conditions 

u velocity 

w wall 

0.5 point in a jet flow where the velocity is ^^max ■*" ^min) concentration is 

2(pi,max ^i,min) 

A bar over a symbol denotes a time mean value, and a prime denotes a fluctuating 

value. 

A double prime denotes evaluation at the reference temperature given by equa- 
tion (33). 


THEORY 

The method used in the present investigation to solve the boundary-layer equations 
for compressible turbulent flow is similar to that used by the authors in previous work 
reported in references 36 and 40. An implicit finite-difference method is used to solve 
the nonlinear partial differential equations for the conservation of mass, momentum, and 
total enthalpy for the mean flow. Details of the numerical procedure and the computer 
program for conventional boundary -layer flows are available in reference 41. The mod- 
ifications to the previous computer program required to solve the present problem are 
described and listed in FORTRAN IV in the appendix. 

In this section, the eddy-viscosity and mixing-length expressions developed previ- 
ously by the authors for calculating turbulent boundary layers without slot injection will 
first be reviewed briefly. The extensions and modifications of these expressions to cal- 
culate the rapidly changing profiles of velocity, temperature, and concentration just 
downstream of the slot and in the subsequent relaxation region (where the profiles relax 
to those for an undisturbed boundary layer far downstream of the slot) will then be 
described. In accordance with previous results for compressible turbulent boundary 
layers (refs. 34 to 36, 40, and 42, for example), kinematic scales are used in the eddy- 
viscosity and mixing-length models. Mixing- length relations determined from incom- 
pressible flows are then applied directly to the present problem. 
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Turbulent Flux of Momentum and Enthalpy 

The turbulent-flux terms in the boundary -layer equations are modeled by eddy- 
viscosity and turbulent Prandtl number concepts, where 


Tt 



and 


q-j. - 


e 8h 
Npr,T 


( 1 ) 


For the boundary -layer flows treated previously (refs. 36 and 40), the eddy viscosity was 
formulated according to Prandtl's mixing-length relation 


e = pj?" 
u 


9u 

9y 


( 2 ) 


where generally 


^u 

6 





(3) 


The quantity within the brackets of equation (3) is Van Driest's damping function, which 
has been modified (ref. 36) to account for the effects of mass transfer on the viscous sub- 
layer. For conventional boundary layers, the function f(x,y/6) consists of two parts. 
These two parts apply, respectively, to the near-wall and far-wall regions, where the 
function was assumed to be (ref. 36) 


(4a) 

g 0 . 3 ) (4b) 

where the form for ff used herein is 

ff = 0.265 - 0.196Hin + 0.0438Hin (4c) 



Equation (4a) is Prandtl's wall function, which seems to be universally applicable. Equa- 
tion (4c) accounts partially for the effects of pressure gradient on the outer, or wake, 

portion of a boundary layer. The limiting values of fjj and ff at ^ = 0.1 and 0.3 

are connected by a straight-line segment. In the present calculations, K = 0.4 has 
, been used throughout. 


The basic concept of eddy viscosity, as expressed by equations (1), is known to be 
faulty (see ref. 43, for example), particularly when the velocity profiles have maximum 


11 


or minimum values within the boundary layer, as in slot-injection flows. For these 
flows, equations (1) require that when = 0, then T.p = 0. This condition is not in 

agreement with experimental observations. An analysis in reference 43 of several sets 
of experimental data, including a wall-jet flow, has indicated that the addition of another 
term proportional to the product. 



would impart the correct behavior to the expression for turbulent shear. This additional 
term is not used herein because its effect is generally small and also because values of 

v'^, which are not generally available, would be required. 

Considerable success has already been achieved (refs. 14, 37, and 38) in the cal- 

A V 

culation of turbulent boundary layers with slot injection for — > 20 by using simple 

s 

eddy -viscosity and mixing-length relations similar to those of equations (1) to (4). How- 

Ay 

ever, for applications of interest herein, the initial relaxation region where — < 20 is 

s 

critical, not only regarding predictions for the extent and magnitude of thermal protec- 
tion, but particularly with respect to the validity of skin-friction predictions. In applica- 
tions of slot cooling to hypersonic cruise aircraft, for example, the magnitude of reduc- 
tions in skin friction may determine the overall feasibility of the slot-cooling system. 

The following models have therefore been developed primarily to provide improved pre- 
dictions in this near-slot region. 

Modified Mixing-Length Expressions 

The type of flow to be computed is illustrated by the schematic sketches of figure 1. 
Injection occurs at x = Xq from a tangential slot of height s with the velocity profiles 
specified in the slot and across the initial boundary layer of thickness 6 q. If the slot- 
lip thickness t is appreciable compared with the slot height s (greater than about 
0.1s), the profiles just downstream of the lip would be affected by local separation and 
reverse flow caused by the thick lip. Since the present method is based on the boundary- 
layer equations, the flow in the separated region cannot be computed accurately. Hence, 
the initial profiles used to start the calculation should be located downstream of any sep- 
arated region. For the present purpose, this downstream distance can be taken as about 
2t where the input velocity profiles would be specified, as indicated in figure 1(a). The 
initial boundary-layer flow mixes with the injected jet flow by turbulent diffusion or 
mixing processes. Far downstream, the velocity profiles tend to approach those of an 
undisturbed turbulent boundary layer. The corresponding development of concentration 
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profiles of air with foreign-gas injection at the slot is depicted in figure 1(b). Note that 
with this type of flow, the mixing region and concentration profiles could readily be mea- 
sured by suitable experimental procedures. Presumably, the mixing process is com- 
pleted when the concentration of air at the wall approaches the free-stream value. 

In figure 1(c), possible mixing-length distributions are shown. Three general 
objectives of the present mixing-length expressions are (1) generally accepted distribu- 
tions of mixing length would be used at the initial input station just downstream of the 
slot location; (2) this initial distribution would be modified in a reasonable and realistic 
manner to represent the early development and spread of the mixing region; and (3) the 
calculated relaxation of the jet and mixing regions of the flow to the undisturbed boundary 
layer far downstream would govern the corresponding mixing-length distributions. The 
mixing-length expressions developed to achieve these objectives and the resulting distri- 
butions of mixing length will be described in the following sections. In order to accurately 
define and "track" the mixing region downstream of the slot, the species- conservation 
equation is utilized. 

In the present computer program, the foreign species is assumed to have the same 
molecular weight and viscosity as air. The present calculations for the concentration of 
a foreign species then represent the physical behavior of a trace species mixed with air. 

Conservation of species .- For no chemical reactions, this equation is written for 
species i as 


__ 8C. 8Cj 

pu — =• + pv 

8x 8y 



8y 


(pv)'c: 


(5) 


A turbulent diffusion coefficient and a diffusion mixing length may be defined 

by the relations 


- 

-(pv)'c; = pDrp — = pi^ic 


8u 

8Ci 


9y 


( 6 ) 


With the use of equation (2), the turbulent Schmidt number Ng^ q> is then simply a ratio 
of momentum to diffusion mixing lengths 


NSc,T = 


e 

PDt 


^u 

Ic 


(7) 


In terms of the concentration profile 


g = 



( 8 ) 
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equation (5) can be written as (where Cj^e assumed constant^ 




9g 


9x 9y ayi^Ngj. ■ Ng^^^y9yJ 
The boundary conditions for g used in the present solutions for y = 0 are 


( 9 ) 


S ~ Sw 

or, for zero mass transfer at the surface 



At the outer edge of the concentration field, — Cj hence, for y = yg 

g = 1 


(10a) 


(10b) 


(10c) 


The requirement of equation (10b) is based on the basic relation between diffusional 
velocity of species i and the concentration gradient. Hence, for zero mass transfer at 


9q\ 

the surface, which is the only boundary condition considered herein, 1 1 = 0. The 

\®y/w 

computation procedure for satisfying this boundary condition is described in the appendix. 


Equation system .- Equation (9) is solved as part of the system of equations for the 
conservation of mass, momentum, and total enthalpy by the same implicit numerical 
method used in reference 36 and given in detail in reference 41. These latter equations 
will not be repeated here since they are identical to those shown in references 36 and 41 
provided that the molecular and turbulent Lewis numbers are assumed to be unity. This 
assumption is used herein. For flows where the turbulent flux of momentum and energy 
is much larger than the corresponding molecular flux, the molecular-transfer terms can 
be neglected entirely except in the near-wall region. The assumption of unity for the 
turbulent Lewis number is based on the experimental data of reference 44, where the 
linear relation between Cj^ and © showed that N^g rp ~ 1.0 for air -air and hydrogen- 
air mixing of subsonic coaxial streams. 


The computation procedure of references 36 and 41 must be changed, however, when 
the recovery temperature of the wall is required. The boundary condition at the wall 
Hw(x) for a specified wall-temperature distribution is replaced by the condition 

= 0 (lOd) 



14 



and the energy equation is solved directly in terms of ? = H/He rather than in terms of 
0. In this way, the values of are obtained directly from the numerical solutions. 

The procedure is the same as that used to satisfy the boundary condition (10b) and is 
described and listed in the appendix. 

The equation of state to be used with the new system of equations would, in general, 
account for the mixing of different gases since the instantaneous properties of the gases 
are related by 

( 11 ) 

i 

When the injected and free-stream gases have the same molecular weight as in the pres- 
ent computer program, 

mi = m (12) 

The mean properties are then related by 

P=^(pT+]^) (13) 


The correlation term p^T’ has been neglected in all solutions on the basis of the results 
and discussion in reference 45, which indicates that for stream Mach numbers up to about 
9, this term would decrease p by only a few percent. 

Calculation of initi al m ix ing region .- Since the system of equations to be solved is 
parabolic, initial conditions for all dependent variables are required. These initial val- 
ues for velocity, concentration (illustrated in fig. 1), and enthalpy are specified as profile 
shapes or functions of y at the initial station Xq. That is, the functions F(xojy), 
g(xo,y), and 0(xo,y) (or ?(xo,y) for boundary condition 10(d)) must be specified. If 
possible, these initial profile shapes should be based on experimental data for the particu- 
lar configuration of interest. In the absence of experimental data for the external flow 
just upstream of the slot and the flow within the slot at its exit, suitable approximations 
for these F and 0 profiles can usually be provided on the basis of general knowledge 
of two-dimensional (or axisymmetric) boundary-layer and channel flows. However, it is 
important to include any upstream '’history'' effects (that is, the effects of static pressure 
and wall-temperature gradients) in these initial profiles for F and 0. If these effects 
are thought to be significant, finite -difference solutions for the upstream flow should be 
used unless reliable experimental data are available. 


15 



With the present definition of g (eq. (8)), the initial concentration profile is sup- 
plied as essentially a step function, with g = 0 at the slot exit and g = 1 for the 
external flow. Thus, for example, in a problem with foreign-gas injection from the slot 
into an air boundary layer, the concentration Cj represents the concentration of air in 
the mixture of air and foreign gas. (See fig. 1.) The molecular weight (and density) of 
the mixture is then governed by the general equation of state (eq. (11)) according to the 
local values of Cj determined from the solution of the complete system of equations. 

When the injected and external gases are homogeneous, as in the present computer 
program, the initial profile for g would be the same as that with foreign-gas injection, 
but the conventional equation of state (eq. (13)) is used and the physical interpretation of 
Ci then corresponds to the concentration of a "trace" species which can have no direct 
physical effect on the velocity and temperature field. The "edges" of the mixing layer 
are then determined at any station from the computed values of g by specifying (in the 
present calculations) the inside edge nearest the wall as the value of y where G = 0.01 
and the outside edge as the value of y where G = 0.99. In the same way, the nominal 
center of the mixing region is specified as the value of y where G = 0.5. (See 
fig. 1(b).) 

Mixing length in free turbulent mixing region.- The value of the concentration 
mixing length within the mixing region between the injected flow and the initial 

boundary layer is assumed to be proportional to the height of the concentration mixing 
region 


^c,m - ^m^ (1“^) 

where W = (y)G=0 99 “ ^^^G=0 OT values of deduced directly from experi- 

mental data are available; hence, the assumption is made that 




u 


(15) 


which is evaluated from velocity data for two-dimensional mixing layers that simulate as 
closely as possible the present initial conditions for g, that is, a step function with 

F = 0 and 1.0 on the two sides of a splitter plate. The assumption that ^ is based 

W Ou 

partly on the results of reference 46 (for axisymmetric jets of water), where by the 


method of Tollmien (ref. 47, p. 412), the concentration mixing length was 


I 


^ 0.5 


0.22 


and the velocity mixing length was ^ a 0.21. ^The parameter rQ 5 ^or yQ 5^ 


IS 
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a more reliable measure of the mixing-layer height than the total height of the mixing 
layer for free turbulent mixing with quiescent surroundings because of the asymptotic 
approach to zero velocity for this type of flow.^ With this assumption, the value of the 
mixing length in the free-mixing region between the jet flow and the exterior flow is 
(from eqs. (7) and (14)) 

^u,m = Nsc,TamW (16) 

This procedure then provides calculated values of a velocity mixing length for the present 
slot-flow problems, where the initial velocity profiles are not uniform and, therefore, the 
values of 6^ rn cannot be determined directly from the velocity profiles. 

To account for the effect of different scales in the outer boundary layer and the 
injected jet on the mixing length within the mixing region proper, the value of Zu,m is 
assumed to apply only at the center of the mixing region, denoted by y^ r “ ^G-0 5 
ures 1(b) and 1(c). The method used to obtain the coordinates of the other "pivot" points 
for the mixing-length distributions in the three zones indicated in figure 1(c) will be pre- 
sented in subsequent sections of this report. 

Experimental values of the turbulent Schmidt number have been measured in sub- 
sonic slot-injection flows, as reported in reference 38, where it was concluded that 
Nsc,T = 0.5 ± 0.2 was representative of the majority of data points. Experimental val- 
ues of Nsc,T’ given in reference 46, for axisymmetric air and water jets vary from 0.67 
to 0.80. Also in reference 46, the value of the mixing-length ratio (which is 

Nsc,T i*^ ih® present approach; see eq. (7)j was 0.81 as determined again by the method 

of Tollmien. The data for subsonic coaxial jets in reference 44 gave average values 
for Npj, <p of 0.6 and 0.85 for air-air and hydrogen-air mixing, respectively. As 
mentioned previously, these same data showed that NLe,T which requires that 

^Sc,T “ ^Pr,T- results given in references 14 and 38 for concentration profiles and 

effectiveness as calculated by the finite-difference method of reference 37 were in better 
agreement with data when a linear variation for Nsc,T was used ^rather than constant 

Nsc,T O’ 5 I’O) that varied from 1.75 at the wall to 0.5 at the boundary-layer edge. 
Since this particular linear variation of Ng^. p is not representative of most data 

(ref. 38) , it was speculated in reference 38 that the improved agreement with data 
obtained with this linear variation was caused by compensating errors in the mixing- 
length model. From these brief comments on turbulent Schmidt number, it is apparent 
that little justification exists for using a particular value or functional variation for 
Nsc,T’ Fov the present solutions, constant values of Ng^, q, of 1.0 and 0.8 have been 
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used. The molecular Schmidt number, Ngg has been taken as 0.7. To be consistent 
with the assumption that NLe,T = the relation that Ng^ -p = Np^ P should be used. 

However, in the present solutions, = 0.9 was used throughout partly in an attempt 

to isolate the effect of changes in Ng^ -p on mixing length, and partly because the value 
Npr,T = 0.9 is representative of typical values used in previous calculations. (See 
refs. 34, 35, and 37, for example.) 


Values of a^n used herein are based primarily on data obtained for two- 
dimensional, free turbulent mixing layers. Although these values are found to vary 
considerably, reasonable limits can be established from data on jets and wakes. 


Thus, from the data of reference 48 for a half-jet, the average value of 
^u 

— = 0.15 for the stream side of the jet. The thickness of the low-velocity side of 

yF=0.99 

the half-jet to the F = 0.025 point was about 1.4 times larger than that of the stream 

side, and the average I was smaller; therefore, ~0.06 for this side (the 

^F=0.025 

F = 0.025 point is used to represent the "edge" of the low-velocity side of the jet because 
of the extremely gradual approach to F = 0 on this side of the jet). If the total thickness 
of the half-jet is taken as the sum of the yp_Q 99 and yp_Q Q25 values used above, the 

average value for am is — = 0.05 where 6 is the thickness of the jet between the 

6 

points where F = 0.025 and 0.99. 


For the plane symmetric jet (ref. 49), the value of 

the center line and edge are excluded. For these data, the value of 

therefore, =0.12 or - — =0.06. 

yF=0.025 ^F=0.025 



0.25 if values near 
^F=0.025 


yo.5 


2 . 1 ; 


For the two-dimensional wake-flow data of Townsend (see ref. 47, p. 394), an aver- 


age value of velocity mixing length (from Tollmien's method) was 




yo.5 


0.4. Therefore, 


— =0.1, since 6 = 2yp.^^ 


F=0.99 y 


yo.5 


F=0.99 


0.5. This value of is in agreement 


with Schlichting's data for wake flow (ref. 50, p. 692) where 
method. 


y =0.09 


by the same 


The values used for a^n in the present solutions were generally in the range of 
0.05 to 0.12. This range of values is based on the experimental results discussed 
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previously; that is, the upper limit is obtained from the symmetric jet data of refer- 

ill / 

ence 49 by assuming that applies directly to in the present slot flows, 

yF=0.025 ' 

and the lower limit is based on the value for the total width of the half jet (ref. 48). 

Mbdng-length distributio n at slot exit and across initial boundary layer .- The 
mixing-length distribution at the slot exit is based on established values for fully devel- 
oped channel flows. In reference 51, for example, it was shown that a mixing-length dis- 
tribution of the form 


I = 




(17) 


yielded accurate velocity profiles and friction-factor equations for fully developed pipe 
and two-dimensional channel flows ^for pipe flows, s is replaced by the pipe radius, 

r^y For channel flows, values of k increased from 0.13 to 0.16 as the Reynolds num- 
ber, based on the average velocity and the hydraulic diameter decreased from about 
9 X 10^ to 2 X 10^. Also, equation (17) gives PrandtPs wall relation (eq. (4a)) for small 
y, with K increasing from about 0.35 to 0.43 for the variation in k noted previously. 


For the present solutions, equation (17) is approximated by a simple two-step ramp 
function. For the near-wall region, the Prandtl slope is used (eq. (4a)) and in the center 
region of the jet, is assumed constant, as given by 



s 




(18) 


where for the present solutions ai was assumed to be 0.14. This value is obtained 

J 1 

from equation (17) evaluated at y = — s with k = 0.14. When the slot flow is fully 

2 

developed turbulent channel flow, this factor could be varied with channel Reynolds num- 
ber based on the results of reference 51. On the other hand, if the slot, flow is believed 
to be laminar, the value of aj would presumably be reduced considerably. 

The mixing-length distribution across the initial boundary layer is also a two-step 
ramp function with equations (4a) and (4c) used in the near- and far-wall regions, respec- 
tively. However, the limiting values of y/6 that are used for these relations in conven- 
tional boundary-layer calculations are not applied herein at the initial input station. 
Instead, the limiting value of y/6 is determined from the intersection of the straight 
lines given by equations (4a) and (4b). 

The coordinates of the pivot points at the initial input station x = Xq are then 
given as follows (see fig. 1(c)): 
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Point 


y 


i 


© 

2K^j 



© 


f 

^ (19) 

© 

^G=0.5 

Nsc,Tam(yc,f - yc,n)^^ 

© 

, ffSo 
s+t+ ^ 

ff«o 

J 


Possible alternative coordinates for point that are more closely related to the slot 

geometry are 




y = s + — t 
2 

and 

= 0.1t 



While these coordinates for point are not included in the present computer program, 

they can be utilized by appropriate adjustment of the input g profile. Note that the 

relations for points and(^ account, in an approximate fashion, for a finite thickness 

t of the slot lip. The mixing-length distributions between the coordinate points given 
previously and in the following sections are taken as straight-line segments. It should 
also be noted that throughout the present method, the Van Driest wall-damping function 
(see eq. (3)) is used only as a modifying factor applied to the basic mixing-length rela- 
tions. Since this damping function depends exponentially on the distance y from the 
wall and this function would have little effect on the final mixing lengths used in the 
program except for small values of y or when the density level is small. 


Mixing-length relations for downstream flow field .- The entire slot-injection flow 
field from the initial development of the mixing region to the final relaxation of the flow 
to an undisturbed boundary layer is calculated by utilizing somewhat different mixing- 
length relations in three distinct streamwise zones of the flow that are indicated in fig- 
ure 1(c). The relations used for mixing length in these three zones are assumed to 
depend on the relative values of Zu,m> ^j? ^b* 

The mixing length in the outer region of the boundary layer is intended to pro- 
vide for a continuous adjustment of the initial input value of If (determined from eq. (4c)) 
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as the initial boundary layer mixes with the jet flow. The assumed scale for this purpose 
is 6 - y^. Hence 

- yc,n) (2«) 

where y^ is taken as yQ^ for the present solutions and where ff g is the same 
function used in equation (4b) except that the lower limit of the integrals for and 

0in is yQ_o 01 I’s^ther than y = 0. 

With the basic relations for the mixing lengths ^b Si^en, respec- 

tively, by equations (16), (18), and (20) (which are to be applied, respectively, to the 
mixing region, the jet flow, and the outer boundary layer), the next task is to specify the 
distribution of mixing lengths across the entire boundary layer and the streamwise devel- 
opment of this distribution. The streamwise development of the distribution in zones I 
and II (see fig. 1(c)) is controlled by the magnitude of relative to the neighboring 

values of Zj and Zj^. In zone III an alternate criterion based on the computed value of 
^ is also used. With the Prandtl mixing-length slope K always applied at the wall 
and with Zu,m applied at the center of the mixing region (y at G = 0.5), the y coor- 
dinates of the "pivot" points connecting the three values can then be specified by adopting 
appropriate scaling factors. The relative values of the mixing lengths are used as cri- 
teria to determine the extent of the three streamwise zones. These criteria and the 
coordinates of the corresponding pivot points, identified in figure 1(c), are presented in 
the next three subsections. 

Zone I (initial mixing region): This zone is defined by the following inequality: 


Zj — fu,m ^ ^b 

The y and Z coordinates of the pivot points are then as follows: 


Point 

y 

Z 

© 



© 


} 

© 

^0=0. 5 

Nsc,T^m(yc,f - yc,n) 

© 

* ^f,s^(^ 6 

K jsr 

^f,s(^ - yG=0.0l) 


( 21 ) 


( 22 ) 
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where 6r = 6 q + t + s. The y coordinates of pivot points and include the fac- 
tor s/Sr in order to maintain the ratio of y/6 the same as the initial values for these 
points from relations (19). The use of this factor is suggested by standard procedure in 
boundary-layer calculations where l/b is often assumed to be a function only of y/6. 
(See refs. 34 and 37, for example.) 

Zone II (intermediate mixing region): This zone is defined by the inequality 


7. < ; <7, 

- ‘^Ujm - ‘•b 

The coordinates of the pivot points are then 


Point 

y 

© 

^Sc,T^m/ N 

K Xc,i - yc,n) 

© 

yo=o.5 

© 

[ ^ 6 

s + t + — 

V K /6r 

0 

/I ^ %,s^o\® " yG=0.5 

K / 6o ^^G=0.5 


(23) 


^Sc,T^m^yc,f “ yc,nj 
Nsc,Tam(yc,f “ yc,n) 

%,s(^ ■ yG=o.oi 


> (24) 


k,s(^ ~ yG=0.01 


where point is a possible alternate form for point . It can be seen that in this 

intermediate mixing region, the reference mixing length is still applied at the 

center of the free-mixing region, but also is extended toward the wall to intersect with 
the Prandtl slope and thereby account for the increasing influence of the entire boundary 

layer on the wall region. Also, the alternate form for y at point (4^ could be used 
when s « 6 q and y at point could then exceed y at point in zone II. 

Zone in (approach and relaxation to final equilibrium boundary layer): The crite- 
ria that control the change to the final mixing-length relations used in this zone are based 
on either the final relaxation of mixing-length relations, as given by the inequality 


— ^u,m — ^b 


(25) 
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or by the computed wall concentration according to the relation 


Ca,w^0.85 


( 26 ) 


This latter inequality is assumed to apply during the final stages of mixing, and the 0.85 
factor was used herein to represent the beginning of this condition. Tests are applied 
continuously toward the end of the solution and in the present computer program (see 
appendix), whichever of these relations (25) or (26) is satisfied first is used as the crite- 
rion to start applying the following coordinates: 


Point 


© 

© 


y 

0.16 

0.36 



(27) 


Thus, in this last zone, essentially the same relations are used that would be applied to a 
conventional undisturbed boundary layer. (See eqs. (4).) The slight difference between 
the limits used for g and is of little consequence in this zone. Since this 

zone should account for the final adjustment of the boundary layer to its undisturbed 


"state," the values of C 


A,w 


should approach free-stream values somewhere in zone III. 


Depending on the choices of am and Ngg^'p, criterion (25) may or may not be satisfied 
before Ca,w -* l-O. 


APPLICATION OF METHOD AND DISCUSSION OF RESULTS 


In this section, theoretical predictions for velocity profiles, heat transfer, and 
effectiveness will be compared with experimental data downstream of tangential injection 
from references 21 to 24 and 29. Also, theoretical predictions of skin friction will be 
compared with data obtained by Aubrey M. Cary, Jr., at the Langley Research Center 
with the same apparatus as described in reference 29. 

The slot configurations used in the investigations of references 21, 22, and 29 con- 
sisted of convergent or constant-area channels, so that essentially sonic velocities were 
obtained at the exit of the slots. The velocity at the slot exit was supersonic in the 
investigation of reference 24. Experimental data were generally available at the slot 
exit or slightly downstream of the exit, and these data were used as the initial inputs for 
the solutions. 

Comparisons With Velocity Profile Data at Mach 3.0 for Adiabatic Injection 

The first two test cases to be considered are taken from data of reference 21 and 
computed profiles of velocity, concentration, and mixing-length distribution will be shown 
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in some detail for these cases. The free-stream Mach number was 3.0, and the slot 
walls were inclined 15*^ to the surface of a flat-plate model. The value of t was taken 
as 0.0025 cm (0.001 inch). The temperatures of the free-stream and jet flows and of 
the surface were such that the flow was essentially adiabatic, and velocity-profile data 

were obtained at seven stations from the slot exit to ® 235. Since the only expert- 

O 

mental data available for these flows are the velocity profiles and skin-friction coeffi- 
cients, the theoretical results for concentration and mixing length are included mainly to 
illustrate how these parameters influence the computations for realistic flow conditions. 


Large slot height, — -0.95 

§o. 

.- This test case is the rearward-inclined step slot 

flow reported in reference 21 with s = 1.7 mm (0.068 inch). Theoretical calculations 

and experimental results for F and 6 for the seven stations are shown in figure 2. 

Also shown are the computed distributions for G, I, and y„ f (G = 0.99), y^ „ and 

G=0.5 


^c n ~ O Ol)- Results from three solutions are presented to indicate the effect of 
changes in the mixing-length parameter ajj^ and in the turbulent Schmidt number Ng^ 
The criterion used in these solutions to control the changeover to the equilibrium mixing 
lengths of relations (27) was that of inequality (25) only. 


The input velocity at the slot lip (see fig. 2(a)) is not zero because the input profile 
was assumed to apply at some small distance downstream of Xq. Thus, while the input 
profile was based on available data, as indicated, a finite value of F = 0.3 was used at 
the slot lip (y = 1.7 mm (0.068 inch)) rather than zero which would be required to satisfy 
the no-slip boundary condition on the slot lip itself. Even with this apparent increase in 
velocities in the vicinity of the slot lip, the integrated mass flow from y = 0 to y = s 
was only 0.0154 kg/sec (0.034 Ib/sec) which is about 10 percent smaller than the quoted 
experimental value of 0.0172 kg/sec (0.038 Ib/sec). This discrepancy in theoretical and 
experimental injected mass-flow rates is probably due to some excess in static pressure 
at the slot exit compared with the local free-stream static pressure, as may be inferred 
from data given in reference 21. 

Comparison of the predicted and experimental velocity profiles in figures 2(a) and 
2(d) shows that solution 2 with the smallest value of aj^ gives the smallest values of 
velocity which are in the best agreement with data for x > 11.05 cm (4.35 inches). The 
smaller velocities in solution 2 resulted from the smaller value of am, which caused 
smaller values of mixing length, as shown in figures 2(c) and 2(f). In fact, for this solu- 
tion the mixing-length relations were such that Zu,m < throughout, therefore, zone III 
(see fig. 1) was never attained in this solution since criterion (26) was not used. At 
X = 10.05 cm (3.95 inches), solution 1 gave the best agreement with the data, whereas at 
X = 8.96 cm (3.53 inches) solution 3 gave better agreement. In view of the difficulties 
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involved in obtaining accurate profile data in a boundary layer 2.5 to 5.1 mm thick (0.1 to 
0.2 inch), all three solutions are considered to be in satisfactory agreement with the data. 


The G profiles and edges of the mixing region are shown in figures 2(b) and 2(e). 
Comparison of results from the three solutions shows that solution 2 ^with am = 0.06 
and Ngc T = 1.0^ gave the smallest values of G at the first two stations downstream of 
the lip. However, at the last three stations, this solution gave larger values of G 
apparently because of the significant reductions in both F and G boundary-layer 
thicknesses for x > 12.7 cm (5.0 inches). These smaller boundary-layer thicknesses 
are again caused by the smaller values of mixing length for solution 2 . 

It is of interest to point out similarities between the computed shape and distribu- 
tion of the concentration mixing region, as given by the values of y^ ^ and y^ in 
figure 2 with the corresponding edges of the mixing region from velocity data of refer- 
ence 9. The experimental coordinates of the upper and lower edges of the mixing region 
resulting from tangential slot injection into a flow with a thin initial boundary layer were 

shown in this reference. For — < 0.36 and — < 10, the inclination or angle of the 

U 0 s 

lower edge was considerably greater than that of the upper edge. This same trend is 
apparent in the theoretical results shown in figure 2(b). Also the location of an effective 
origin of the velocity mixing region (determined by upstream extrapolation of the linear 
part of the edges) in the data of reference 9 was located slightly upstream and above the 
slot rather than at the slot lip itself. This same general shift in the effective origin of 
the mixing region is also apparent from the results for y^, j and y^. shown in fig- 
ure 2 (b). 


The locations of the three mixing-length zones defined in figure 1 are indicated in 
figure 2(c) for these solutions. It is seen that the values of a^ have a marked effect 
on the location of the changeover points from one zone to the next. Thus, while the best 
agreement with the downstream velocity profiles was obtained with solution 2 , the failure 
of this solution to reach zone III or final equilibrium mixing-length distributions even at 

the last station = 23 5^ is probably an indication that either am was too small or the 

alternate criterion (26) should have been used. In this connection, it is of interest to 
compare the values of from the three solutions at x = 20.3 and 49.0 cm (8.0 and 

19.0 inches) given in the following table: 


Solution 

^m 

^Sc,T 

Ca ai X in cm (in.) 

equal — 

20.3 (8.0) 

49.0 (19.3) 

1 

0.09 

0.8 

0.81 

0.93 

2 

.06 

1.0 

CO 

.90 

3 

.10 

1.0 

.78 

.92 


25 



Again, ^ values for solution 2 are the smallest, whereas the values for the other 
two solutions approach unity more rapidly; this trend indicates that the predicted mixing 
process is faster for solutions 1 and 3. Thus, the = 0.85 criterion of inequal- 

ity (26) was attained first in solution 1. This result is believed to be the most realistic 
because the experimental velocity profile shapes were nearly the same as those for the 
undisturbed flat plate (see comparisons in ref. 21) by, at least, the x = 36 cm 
(14.0 inches) station. Hence, the values of am and Ngc,T used in solution 1 have 
been selected as reasonable compromise values and will be used in all remaining solu- 
tions, partly on the basis of the behavior and partly because of the good agree- 

ment with velocity profiles from this solution and the data for the first two stations in 

the critical region of — < 10. Also, the changeover criterion given by inequality (26) 

s 

will generally be used in the remaining solutions. 

Small slot height, == 0.20 

.- This test case was also taken from reference 21 but 

for the rearward step slot with s =0.38 mm (0.015 inch). The predicted F and G 
profiles and the mixing-length distributions are shown in figure 3. The mixing-length 
ratios used in this solution were aj = 0.14 and aj^^ = 0.09 and the turbulent Schmidt 
number was Ng^ = 0.8. The experimental velocity-profile data are also shown in the 
upper part of the figure. 

Since the quoted value of injected mass flow was 0.0172 kg/sec (0.038 Ib/sec), the 
same as for the large-slot-height test with s = 1.7 mm (0.068 inch), the experimental 
injection pressure would have to be increased by the ratio of s values or by about 4.5. 
Hence, the velocities just downstream of the lip would be considerably larger than sonic 
because of the sudden expansion of the slot flow. The input -velocity profile is again 
based on experimental values which do indicate (fig. 3(a)) that a large and sudden expan- 
sion apparently took place just aft of the lip as evidenced by the outward displacement of 
the jet flow and the large ratio of jet to external velocities of about 0.8 at the first station. 
Since the external Mach number was 3.0, the injected jet flow was supersonic at about 
Mach 2.4 at this first station. Since the injection channel was constant area or conver- 
gent, the pitot probe would therefore have been slightly downstream of the lip. With the 
present theoretical restriction of constant static pressure, the integrated theoretical 
injected mass flow was only 0.0099 kg/sec (0.022 Ib/sec) or about 40 percent lower than 
the experimental value. 

In spite of this discrepancy in injected mass-flow rates, the computed velocity pro- 
files are considered to be in reasonable agreement with the data. This acceptable agree- 
ment indicates that reasonable predictions of velocity distributions can be obtained even 
when the static pressure and angle of the injected flow deviate considerably from the 


26 



theoretical assumptions involved in the boundary-layer equations ; providing initial pro- 
files are available. 

The mixing-length distributions shown in figure 3 indicate that zone III (and the 
associated final equilibrium mixing lengths) was reached upstream of the station where 
X = 20.4 cm (8.05 inches) in this solution. Actually, the criterion of ^ - 0.85 was 
reached at x = 20.3 cm (8.00 inches), just slightly upstream of the station where pro- 

A V Ax 

files are shown. This corresponds to = 310 and ~ 40, either of which seems a 

physically reasonable distance to begin the approach to an equilibrium boundary layer. 

For the conditions of this solution, the mixing-length criterion (25) was never reached 

= 1060 or — ~ 90 j, which is probably an exces- 
sive distance to reach equilibrium. Hence, for situations like this one, where s « 6 q, 
criterion (26) is thought to provide more realistic results. 

Skin-friction predictions and comparisons with expe rimental data .- The calculated 
distributions of skin friction for s = 0.38 mm and 1.7 mm (0.015 and 0.068 inch) are 
compared with experimental data from reference 21 in figure 4. At x = 50.8 cm 
(20 inches), data were obtained with a skin- friction balance on the solid plate (no slot) 
and for both rearward-inclined step slots of s = 0.38 and 1.7 mm (0.015 and 0.068 inch). 
These data showed that the 1.7-mm (0.068-inch) slot caused a 5-percent increase in skin 
friction while the 0.38-mm (0.015-inch) slot caused a 3-percent decrease. The other 
data for zero injection are based on measured values of 6. Comparison of the theoreti- 
cal predictions with these no-injection data indicate that injection reduces the skin fric- 
tion over the entire plate for both slot heights. The greatest reductions were predicted 
for the 1.7-mm (0.068-inch) slot. Significant reductions were also predicted for the 

^ y 

smaller slot, but only up to x “ 10. 16 cm (4 inches) or — ~ 40. The magnitude of 

s 

actual reductions in Cf that may be realized for mismatched pressure conditions like 
those for the smaller slot should probably be determined from experimental data since 
the expansion or compression wave systems would probably modify the present velocity 
profiles considerably. 


Comparisons With Heat -Transfer Data at Mach 6 for 



0.66 


Sonic injection velocities .- The next test case to be considered is based on data 
from references 22 and 23. The model was an axisymmetric center body mounted on 
the center line of an axisymmetric Mach 6 contoured nozzle. The diameter of the model 
in the test region was 11.74 cm (4.623 inches) and the initial boundary-layer thickness 
just ahead of the slot was about 2.03 cm (0.8 inch). The data to be used here were 
obtained with a rearward-facing tangential slot of 0.25-mm (0.01-inch) height. The area 
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distribution of the approach channel to the injection slot was convergent, so the injection 
velocity was sonic. 

The mass-injection ratio X was varied from about 0.06 to values greater than 
unity by increasing the injection pressure. The mass-injection ratio is defined herein 
as 



In terms of an average mass ratio, 


Xav - 


\ ^e'^e /av 


where Xav “ the ratio of the static pressure at the slot exit to the local-stream static 
pressure is 


Pj,o 


1 + 


y - 1 ^>.2 


1/2 


m: 


e T 


t, 3 ,o 


Pe ^j ,0 I y ~ ^ m2 ^ 


'^av 


(29) 


3,0 


'^t,j,o ^ 


Then with the conditions of references 22 and 23 I Mj = 1.0, Mg - 6, and - ~ 0,661, 


^3,0 

Pe 


~ 12.7Xav 


(30) 


and in order to satisfy the present limitation of constant static pressure the value of Xav 
should be about 0.08. 

The input profiles used in the solutions are shown in figure 5. These profiles are 
based on data given in reference 23 for y ^ 0.25 mm (0.01 inch). The profiles in the 
slot region were adjusted to give X = 0.068 so as to satisfy approximately the require- 
ment of ~ 1.0 from equation (30). All data reported in references 22 and 23 were 

Pe 

T 

obtained with T^, approximately constant at room temperature. Hence, — — = 0.66 

^t,e 


represents the mean of the experimental range in Tj.^g of 400° K (720° R) to 500° K 

(900° R). A solution designated in the following discussion as the "heat-transfer solu- 

T 

tion” was then obtained with — — = 0,66. 

Tt.e 
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The predicted heat-transfer distribution from this solution is plotted in figure 6 in 
terms of the correlating parameters from reference 22 |l - ^ 

The reference heating-rate values ^ were computed from the same reference- 
enthalpy expression used in reference 22: 



where 


•aw 


= w Tt,e - 


+ T. 


(32) 


and o) = 0.89 was used herein to obtain q,j„ The gas-property values in equation (31) 

W,I 

are evaluated at the reference temperature (ref. 22) 


T ' = 0.5T 


w 


1/3/ l/3\ 

+ 0.22Npr (0.5 - 0.22Np'^^ lTt,e 


(33) 


It was stated in reference 22 that the measured heating rates for zero injection were within 
10 percent of the values from equation (31), so the computed values of q^ j. are con- 
sistent with the experimental values. 

Also shown in figure 6 are the experimental data (shown as the hatched band) and 
the straight line which correlated the data in reference 22 for X > 200. The predicted 
heat-transfer distribution is within the spread of the data for 150 < X < 600. The pre- 
dicted results are greater than unity for X < 150 because the input wall temperature of 
Tw = 294° K (530° R) was larger than the effective recovery temperature in this region 
and therefore the calculated heating was negative (heat transfer from the wall to the flow). 
The disagreement between predictions and data for X > 600 is probably within the 
experimental uncertainty band of the data. At X = 800, for example, the theory predicts 
a heating rate that is only 4 percent larger than the largest measured value. (Some of 
the spread in the data may be caused by different values of X; however, the values of X 
for the original data points in reference 22 were not identified.) 

A parameter of particular interest in the application of slot-cooling techniques is 
the effectiveness, which is defined in reference 22 as 


V = 


'^t,e ~ ^aw 
T - ^ 

t,e t,j,o 


(34) 
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It can be seen that rj provides a measure of the thermal protection afforded by injection 
since values of rj near 1 indicate that the local recovery temperature at the wall is still 
close to the jet total temperature. However, since the recovery temperature for the jet 
flow may be less than the jet total temperature, the values of rj from equation (34) could 
exceed unity. 


Values of T^w hence tj are computed in the present method by replacing the 
boundary condition Hw(x) (see ref. 41) by (— ) = 0. The value of required to 

\^y/w 

satisfy the new boundary condition is then computed directly from the differential equation 
for K. Since the surface heat transfer is zero when (42) = o, the resulting values of 




^aw 

Tt,e 


for constant Cp. Details of this procedure, which is also used to determine 


^A,w boundary condition ( 10 b), are presented in the appendix. 


Another solution, designated the "adiabatic-wall solution,” was obtained with this 
procedure for computing Taw with all other input and boundary conditions identical 
to those of the heat-transfer solution. The resulting values of 77 are plotted against the 
correlating parameter X in figure 7. Shown for comparison in the figure are the exper- 
imental data and straight-line correlation from reference 22. The computed values of 77 
from the adiabatic-wall solution are considerably larger than the experimental values of 
77 for X > 200; hence, the computed recovery temperatures are smaller than those used 
in reference 22. As mentioned previously, the values of T^^ used in reference 22 to 
determine 77 were calculated by utilizing measured values of q^ in equation (31). 

When this same procedure for computing 77 is applied to predicted values of q^ from 


the heat-transfer solution 


^obtained with 


Tw 

Tt,e 


0.66 


the 77 


distribution, shown as the 


dashed line in figure 7, is obtained. As would be expected, this distribution of 77 is in 
agreement with the data since the predictions for q^ agreed with the data. 

The values of Taw^T^ g computed directly in the adiabatic-wall solution and com- 
puted from equation (31) with q^ from the heat-transfer solution are shown in figure 8 . 
The latter procedure yields values of Taw/Tt,e most 13 percent larger than 

values computed directly in the former solution. The maximum difference in temperature 
is about 44° K (80° R), which can be easily measured in a facility with long test times. 

The larger values of Taw/^t^e questionable because even for a simple flat-plate 
flow, it is known that equation (31) is not correct when dT^/dx is sufficiently large. 

(See refs. 31 and 32.) In the present situation, the wall temperature increases by about 
83° K (150° R) in a distance of 2.5 to 5.1 cm (1 to 2 inches). Also, the injection and 
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mixing of a cold layer of air near the surface would cast further doubt on the use of equa- 
tion (31) to compute Taw Of course, direct measurements of are required to 

verify the theoretical predictions of the adiabatic-wall solution, which indicate that 
may not reach the nominal equilibrium values until Ax/s exceeds several thousand. 

Such measurements are available (ref. 29) and will be compared with predictions in a 
subsequent section. 

As a matter of interest, the computed values of from the heat-transfer solu- 

tion are shown in figure 8. Again the final equilibrium value for Ca,w unity may not 

A •V' 

be reached until — > 2000 to 3000. 
s 

Supersonic injection velocities .- Surface heat -transfer measurements are reported 
in reference 24 for supersonic injection from. a tangential slot into a 2.54-cm-thick 
(1-inch) turbulent boundary layer. The free-stream test conditions and the model con- 
figuration were the same as those of reference 22 except the slot height s was 5.6 mm 
(0.22 inch) and the injection velocity was supersonic at Mj == 2.3. Accordingly, the the- 
oretical solution for these conditions was obtained with the same initial profiles for F 
and ^ shown in figure 5 except the maximum velocity ratio at the jet exit was increased 

/ Tt 

to 0.617 corresponding to values quoted in ref. 24 of = 2.3, = 0.66, and 

V '^t,j,0 

Mg = 6.2^ and s was increased to 5.6 mm (0.22 inch). 

The predicted heat-transfer distribution in the form of 1 - is plotted against 

^w,r 

X in figure 9. The value of ^ was again computed from equation (31) with the same 
stream conditions and co = 0.89 as used for the previous case. The experimental data 
from reference 24 are also shown in the figure as the hatched area. Data points for 

from 0.139 to 0.248 are included. The theoretical prediction is near the upper edge of the 
hatched data band. However, all data points, except one, in the range of of 0.139 to 

0.248 are near the upper edge of the shaded band in figure 9. Since A = 0.2 was used in 
the solution, the agreement between predicted results and data is reasonable. 

Comparisons With Recovery Temperature and Skin- Friction 

'^t i 

Data at Mach 6 for — — ~ 0.65 

Tt,e 

The data considered in this section are from reference 29 where direct measure- 
ments of Tg^^ downstream of a two-dimensional tangential slot were reported. The 
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free-stream Mach number was 6.0, the free-stream and jet total temperatures were 
approximately 472° K (850° R) and 289° K (520° R), respectively, the exit velocity at the 
slot was sonic, and the initial turbulent boundary layer just upstream of the slot was about 
5.1 cm (2 inches) thick. The ratio of measured slot mass-flow rate to calculated free- 
stream mass-flow rate was varied from 0.06 to 1.6, corresponding to pressure ratios of 

— -0.8 to 20 from equation (29). The model was mounted on the floor of the Langley 
Pe 

20-inch Mach 6 tunnel. 


The initial profiles used for the solution are shown in figure 10. The F and ? 
profiles are based on both the experimental data obtained on the wind-tunnel floor without 
the slot present and on the indicated experimental values of Tw,o/Tt,e» and t. While 

the profiles at the slot exit were not measured, they are based on known properties of 
channel flow. The shape and level of the ^ ‘ profile at the slot exit represents an attempt 
to satisfy measured boundary conditions on the upper surface of the channel and a short 
distance upstream of the slot exit. Also, the value of Tj- j oy/Tt g used to calculate 77 
was 0.65 rather than about 0.61 as used to reduce the data in reference 29. This larger 
value of Ttj^oy^Tt^e was used herein because of increased values of Ttj^o measured 
recently by Cary in the jet exit. These larger values of Tt were apparently caused 
by heat conduction through the steel lip of the slot. The integrated slot mass -flow ratio 
for the profiles shown was X = 0.065 which corresponds approximately to a matched 


pressure condition where 



Pe 


in accordance with this limitation of boundary-layer 


theory. 


The computed distribution of effectiveness 77, defined by equation (34), is plotted in 
figure 11 against the same correlating parameter X used previously. The experimental 
data from reference 29 are shown as the hatched band. Agreement between the predic- 
tion and data is good over the entire range of the measurements up to X 800, although 
the prediction is somewhat higher than the data for 100 < X < 200. The generally good 
correlation of the data with the parameter X for the large range of experimental mass 
ratios from X = 0.06 to 1.6 indicates that the corresponding range of jet-to-stream 


Fi 

pressure ratios of — - 1.0 to 20 has little effect on recovery temperatures. From the 
Pe 

good agreement between the theoretical prediction and the data, it follows that for this 
range of conditions, the surface temperatures are determined mainly by the turbulent 
mixing phenomena as predicted by the present boundary-layer theory. 


The correlation from reference 22 is also shown for comparison in figure 11. At 
large values of the correlating parameter X, the effectiveness from the direct measure- 
ments of reference 29 and as predicted by the present theory is as much as 100 percent 
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larger than the correlation of reference 22. These results show that Taw cannot be 
computed simply from measured heating rates and flat-plate heat-transfer relations. 

Direct measurements of sMn friction have also been obtained recently by Cary at 
three stations downstream of the slot on the same model used for the recovery- 
temperature data of reference 29. Skin-friction balances similar to those described in 
reference 52 were used to obtain these new data. Predicted values of Cf from the solu- 
tion are compared with the data in figure 12. The agreement is good, in spite of the large 
lip thickness (t/s = 1/3) used in the experiment. Two more or less compensating effects 
may be present in this experiment. These are a possible low turbulence level in the slot 
flow caused by the very short approach configuration with its large accelerations (see 
ref. 53) and the thick lip which presumably would tend to increase turbulence levels down- 
stream of injection. These effects can be roughly accounted for in the present theory by 
adjustment of the initial velocity profile in the slot and by using different values of aj. 
Such adjustments are not warranted without data that would be required to incorporate 
properly these effects of turbulence level on the mixing-length models. 

Also included in figure 12 are measured values of skin friction on the tunnel wall 
without the slot. These measurements were obtained with the same instruments, at the 
same x locations, and for the same free-stream conditions as the slot-injection data. 
Comparisons between these data and the data with injection show that the skin friction is 
reduced significantly by slot injection under these conditions up to Ax/s values of at 
least 70. An average reduction of approximately 60 percent is attained in the near-slot 
Ax 

region of < 20. It is also concluded that reliable predictions of level and trends in 
s 

skin friction are obtained from the finite -difference solutions with the present mixing- 
length equations and constants. 


CONCLUDING REMARKS 

An implicit finite-difference method has been used to solve the partial differential 
equations for compressible turbulent boundary layers with tangential slot injection. The 
turbulent-flxix terms have been modeled with eddy-diffusivity and mixing-length concepts. 
The species-conservation equation is used to calculate the concentration field, which, in 
turn, is utilized to provide appropriate scales and criteria for different mixing-length 
models applied in three different zones downstream of the injection slot. When the pri- 
mary and secondary flows are homogeneous, as in the example problems treated herein, 
the species conservation equation governs the spread of a ^^trace” species which has no 
direct effect on the velocity and enthalpy fields, but which serves to define the center and 
edges of the free-mixing region. The mixing lengths for this species concentration field 
are related to velocity mixing lengths by the turbulent Schmidt number. 
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Predictions for velocity profiles, heat transfer, effectiveness, and skin friction have 
been compared with experimental data from four different investigations at stream Mach 
numbers of 3.0 and 6.0 with sonic and supersonic injection velocities. These compari- 
sons show that the turbulent -flux models developed herein provide realistic predictions 
for the entire flow from the near-slot region through the final relaxation region even when 
the boundary-layer limitation of small normal-pressure gradients is violated. This limi- 
tation was violated in various degrees in most of the experimental data considered herein 
because of such factors as mismatched pressures at the jet exit, finite injection angle, 
and thick slot lip. 

An indirect procedure used by previous investigators for calculating the recovery 
temperatures from measured heating rates and flat-plate-correlation equations for heat- 
transfer coefficients underestimates the thermal effectiveness of tangential slot injection 
in hypersonic flows by significant amounts. The flat-plate-correlation equations cannot 
be used in this way because the relaxation process to undisturbed equilibrium conditions 
is slower for recovery temperatures than for heat -transfer rates. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., March 2, 1971. 
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APPENDIX 


MODIFICATIONS TO NONSIMILAR BOUNDARY- LAYER PROGRAM FOR 
COMPUTATION OF SLOT -INJECTION FLOWS 

By Barbara A. Hixon and Dennis M. Bushnell 
Langley Research Center 


The compressible turbulent boundary-layer computer program described in refer- 
ence 41 has been modified for application to tangential-slot-injection flows. The purpose 
of this appendix is to describe these modifications so that the program as presented in 
reference 41 can be used to compute relaxing tangential-slot flows at matched exit- 
pressure conditions. The program can be run with either the adiabatic-wall boundary 
condition or a specified wall temperature and wall -temperature gradient. If ZETWTAB 
table) and DZDXTAB ^d^^/dx/dL) (see ref. 41) are given as input, they will be used; 
otherwise the program utilizes the adiabatic-wall boundary condition. The spanwise or 
g-momentum equation of reference 41 has been used herein as a species -conservation 
equation (eq. (9)). Therefore, it was necessary to discontinue treating g as a velocity 
term (as was done in ref. 41). 


For convenience the complete program listing is given. Portions of the program 
which differ from that in reference 41 are indicated on the listing. The reason for a 
modification is generally identified as being one of the following: 


© 

© 


Additional input and output statements 

Changes in eddy -viscosity expression to account for slot and mixing-region 
flows 


Logic used in computation of adiabatic- and impermeable-wall boundary 
conditions 

The extra input ^indicated by is necessary to compute the concentration pro- 

file, the various pivot points associated with the mixing length for the slot and mixing- 
region flows, and the mixing-length predictions. 


The changes due to the eddy-viscosity function ^indicated by \^Jj are necessary to 

compute the relaxing slot flow and are discussed in the text of the present report. 

To apply the correct boundary condition on the concentration profile for an imper- 
meable wall (zero gradient at the wall, eq. 10(b)), the solution must be obtained in terms 
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of the concentration directly. The algorithm for the computation of the g quantity 
(treated herein as species concentration) is, in the notation of reference 41, 

Sn = Gngn+i + §n (Al) 

where and gj^ are the recursion functions. Since n = 1 corresponds to the wall 
value, the impermeable-wall boundary condition requires that 

g2 =Sl 

and therefore from equation (Al) 

gl = Gjgj + gj 

This result can only be correct for arbitrary gj if Gj = 1 and 

are then used in the recursion relations for G and g at n = 2 
from the wall), with the result (see eq. (31) in ref. 41). 

G 2 = and §2 = ^ 

B2 + C2 B2 + C2 

In the notation of reference 41, the A, B, C, and D quantities are coefficients in the 
general difference equation for g. The values for G 2 and g 2 obtained from equa- 
tions (A2) are used to compute the G^ and g^ quantities (ref. 41) and the zero- 
gradient boundary condition at the wall (eq. 10(b)) is then satisfied. 

When the adiabatic-wall temperature is required, the solution is obtained in terms 
of ? instead of 0 so that again the boundary condition of equation 10(d) can be applied 
directly by the same procedure used for the concentration. These changes are indicated 

by(5)- 

The following is a hist of nomenclature added to the program described in refer- 
ence 41; 

Input: 

Program notation 
Al 

A2 
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Report notation Description 

aj Mixing-length ratio for center of initial 

slot flow, recommended value, 0.14 

am Mixing-length ratio for free -mixing 

region, recommended value, 0.09 


gj = 0. These values 
(the first step away 

(A2) 



Program notation 

Report notation 

Description 

SM 

Nsc 

Molecular Schmidt number 

SNT 

Nsc,T 

Turbulent Schmidt number, should be 

approximately same as turbulent Prandtl 
number, recommended value, 0.8 

TDL 

t/L 

Nondimensional slot-lip thickness 

XXLl 

(x/L) 1 

Nondimensional distance at beginning of x 
step during solution 

XXL2 

(x/L)2 

Nondimensional distance at end of x step 
during solution 

YCDL 

s/2L 

Nondimensional distance to center of slot 

YDDO 

s + t ^f 
6o K 

Initial nondimensional distance to pivot 
point number (see eqs. (19)) 

YDEL 

6 

Local boundary-layer thickness (value of y 
where F = 0.995) 

YDELO 

(s + t + 6 q)/L 

Initial nondimensional boundary-layer thick- 
ness (for entire flow) 

0L 

L 

Reference length, generally 1 cm or 1 inch 

Output: 

YDLG 

yc,r/L 

Nondimensional height to G = 0.5 location 

YSIDL 


Nondimensional height to G = 0.01 location 

YS2DL 

yc,f/^ 

Nondimensional height to G = 0.99 location 


A complete listing of the program for slot-injection flows, with changes to the pro- 
gram of reference 41 identified, is now presented. 
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PROGRAM D3630 ( I NPUT . OUTPUT ♦ TAPES * I NPUT * T APES^OUTPUT ) 

DIMENSION DELETA (350 )*FTAB( 100)1 FI ( 350 ) i F A ( 3S0 ) ♦ FA2 ( 350 > f 
1VA(350).GTAB(100) .G1 (350 > * G2 ( 350 ) » THETA I ( 350 > » THET A2 ( 350 ) ♦ 
2THETAA (350 ) *RHOTAB( 100 ) iRHOEROI ( 350) iRHOEROA (350 ) iZETATABC 100 ) » 
3ZETA(350)iXMBARl ( 350 ) • XMC I RC 1 ( 350 ) ♦ XMSTAR 1 ( 350 ) ♦ XMPR I Ml (350)* 
4PHIR(350) *XM0ARA2 (350) . XMC I R A2 < 350 ) *XMSTRA2 (350 ) *XMPRMA2 (350 ) * 
5CAPG(350) fSMLG(350 ) *eTATAB( lOO) *OUOXTAB (75) iDZDXTAB (75) * 

6 XMBAR2 05C ) • XMC I RC2 ( 350 ) *XMSTAR2(350 ) * XMPR I M2 ( 350 ) * 

7XM0ARA ( 350 ) * XMC I RCA ( 350 ) * XMSTARA ( 350 ) ♦ XMPR I MA ( 350 ) * SUMRER ( 350 ) i 
8F2(350 ) *XLPR(30 ) *VWTAB(75 ) iRMUTAB ( 75) *ABTAB(20 ) *FCFTAB(20) • 
9XMMEG ( 350 ) • XMMEF ( 350 ) * RRUUER(75) * XL ( 75 ) * X t T A0 ( 75 ) * 
1ETA(350)*UEOSTAB(7S) ,RTAB(75) *RERSTA8(75) *YL( 100) *ZETWTAB(75) 
COMMON/EPSOMU/EPSDMU ( 350 ) 

COMMON/TABLE 1/RH0ER02( 350 ) 

C0MM0N/TABLE2/YDL ( 350 ) * RHORHOE ( 350 ) 
common/three/numeta *nmaxf 

C OMMON /FEB I 2 /VW A 

COMMON/FEB 1 1 /PRTTAB ( 20 ) * YOOPRT ( 20 ) * NYP 

COMMON/HDCAPHE/HDCAPhE (350 ) *GEE2 

COMMON/MUUSE/ I USEEMU , MPWEMU 

COMMON/DEC 1 8/YCDL 

COMMON/AGA I N/TDL 

COMMON/DEC22/YS 1 DL ♦ YS2DL 

COMMON/ JAN7/X0 

COMMON/MAY 1 1/Al * A2 

COMMON/ JUNE 1 /SM * SNT 

C OMMON / JUL Y2 /Y DE L 0 

COMMON/HIS/HIS 

COMMON/ IFTC/IFTC 

COMMON/CCW/CCW 

COMMON/ZETW/ZETWTAB 

COMMON/ZETAW2/ZETAW2 

COMMON/FBAR/FBARTAB ( 20 ) * I FBLU* YDDFB ( 20 ) *NFBY 
COMMON/ I WLDMP/ I WLDMP 

NAMEL1ST/NAM1/NUMETA,NMAXF.NMAXG*0ELETA*XK,X10*DELXI0*XITEST* 

1 XISTOP*FTAB*ETATAB*VWTAB*EPSLONE *EP3LONW* GTAB*UEOSTAB * 

2 WEDS2HE *PR . ZETWTAB * XNBAR , RERSTAB ♦ CAPRS ♦ RTAB , J , RHOTAB » 

3SHE * ZET AT AB • HSHE . XL « NUMX * YL * NUM Y ♦ ^ 

4XO*OL*DUDXTAB*DZDXTAR*FRO*NSTEPS* 

5AP*BP.CP*ABTAB *FCFTAr»NFCFAB * PRTTAB ♦YDDPRT» NYP ♦ 

6XLPR. IVEG* IN IT, I USEEMU , MPWEMU , YCDL * TDL * A 1 * A2 , SM , SNT , YPELO * 
7FBART AB , I FBLU , YDDFB * NFBY • I WLDMP 
ZETWTAB(2)^0,0 
ZETATAB(2)=0,0 


I 


Removed ITHETA from NAMELIST 






RH0TAB(2)=0«0 
ETATAB(2 )-0«0 
XLPP (2 )*0«0 
XiOsO.O 

NPRINTsi 
T COUNT 3 I 
IFTC^O 

READ(5.NAM1 ) 

WR1TE(6*NAM1 ) 

OPTION FOR INITIAL DFLTA Xl 


IF(INIT.EQ.O)OELXI=DFLXIO 
IF( INIT.EQ* 1 )DELX1=XI0*1 .E-5 
GEE2*1. 4 


TRANSFORM X TO XI 


HRDCPHE=1 .-UEDSTABC 1 > **2-WEDS2HE**2 

RMUTAB< 1 )=SQRT<HRDCPHE/HSHE ) •» ( HSHE + SHF ) / ( HRDCPHE + SHE ) 
1»RERSTAB(1 >*HRDCPHE 
2/HSHE 

FUNCX 1 3RMUTAB ( 1 ) *UEDST AB ( 1 ) *RTAB ( 1 ) *♦ ( 2* J ) 

XI =xo 

DO 22 1*1 .NUMX 
IF( I .EQ.l )GO TO 15 

HRDCPHE*! .-UEOSTABf f j ♦♦2-WEDS2HE**2 

RMUTAB ( I ) =SQRT ( HRDCPHE/HSHE ) ♦ ( HSHE + SHE > / ( HRDCPHE + SHE ) 
1*RERSTAB( I )»HRDCPHE 
2/HSHE 

FUNCX2 = RMUTAB( I )*UEDSTAB( I )»RTAB( I )*♦ (2*J ) 
SUMFXI=SUMFX1+(FUNCX2+FUNCX1 )/2.» (XL ( I )-XL ( I-l > > 
FUNCXl =FUNCX2 
GO TO 95 
15 SUMFXI*0« 

95 RRUUE R{ I )= CAPRS*RMUTAB ( 1 )#UEDSTAB ( I )*RTAB ( 1 )*♦ (2*J ) 
IF(X10,E0*0. )GO TO 25 
GO TO 31 

25 XlOaRRUUERM )»X1 
31 XITABC I )3XIO-fCAPRS*SuMFXI 

DUDXTA0 ( I )3DUDXTAB( 1 )/RRUUER ( I ) 

IF(ZETWTAB(2)«NE.0. m ZPXTABC I )=DZDXTAB( I )/RRUUE R( I ) 

22 CONTINUE 
XI 1*XI0 


(jO 

CD 


Instfnu of GLEl’ = GTAB(r?) 




O 


IF NO ETA TABLEfTRANSFORM Y TO ETA 
IF(ETATAB(2 ) .NE.G. )GO TO 26 

ZETAWl aZETATABM > 4 : 

IF(ZFTWTAB(2 ) •EQ.O. )r,0 TO 3 4 

DO 202 I*UnUMY 

THETA 1 ( I ) = (ZETATAB( I )-ZETWTAB( 1 ) )/( 1 .-ZETWTAB( 1 ) ) 

202 CONTINUE 

3 ZWsZFTATABU^ 

IF (ZETWTA0(2 ) .EQ.D. )7W«0. 

DO 33 Ial*NUMY 

IF(ZETWTAB(2 ) .NE«0. )GO TO 32 

thetai < n*ZETATAB ( n 

32 RHOTAB( I) = ( ( 1 •-ZW)*ThETA 1 ( I ) -f ZW-UEDSTAB ( 1 )**2 ^ 

!*FTAB ( I )»*2-WEOS2HE **2*GTAB( I )**2 )/( I .-UEDSTABC 1 )**2-WE0S2HE 
2 ** 2 ) 

33 CONTINUE ^ 

DO 34 I*1,NUMY^^ 

RHORHOE (1 ) s I . /RHOTAB ( I ) 

34 CONTINUE 

FUNCYl aRHORHOEM )*RERSTAB( 1 ) 

SUMFETAaO.O 

ETATABI 1 )=0,0 

DO 27 I*2*NUMY 

FUNCY2* RHORHOE ( 1 ) *RERSTAB ( 1 ) 

SUMFETA*SUMFETA+(FUNCY2+FUNCYI )/2.# (YL ( I )-YL ( I-l ) ) 

FUNCYl SFUNCY2 

ETATAB( I )=CAPRS*UEDSTAB( 1 )*RTABU )**J/(2*X10)**XNBAR *sumfeta 
27 CONTINUE 
26 CONTINUE 

constants 

ZETAE=1 • 

FTESTa. 99999 
GTEST= #99999 
THTESTa. 99999 
ETA( 1 )s0. 

SUMRER( l)aO. 


COMPUTE DELT 


ZETAW1=ZETWTAB(i) has changed position in program 
Instead of IF(ZETATAB(2) .0^.0. ) GO TO 32 


Instead of ZETWTAB(l) 
Removed 32 COKTIKUE 



FI (1 1»FTAB( 1 ) 

G1 (1 I=GTABn ) 

IF{?ETATAB(2)tNE.O. IZETA ( 1 )=ZETATAB ( 1 > 

IF (RHOTAB (2 ) »NE*0» JRhOEROI ( 1 )=RHOTAB ( I ) 

Ma2 

DO 10 r»2,NUMFTA 

OFLETA ( I ) ( 1 - 1 ) *DELETA ( 1 ) 

COMPUTE ETA S 

interpolate to get F(ETA)S 

ETA( I JsETAI I-l )+0ELETA (1-1) 

IF(F1 ( 1-1 ).GT..9)M=:I 

CALL FTLUPIETA (I ) «F1 (I ) •MtNUMY *ETATAB , FTAB ) ^ 

W=I S 

CALL FTLUPIETA ( I ) *G1 (1 ) ♦ M ,NUMY* ETAT AB « GTAB ) 

M*2 ^ 

IF(ZETAM-1 )*GT*«9>Ms1 ^ ^ 

CALL FTLUP(ETA< I ) <ZEtA( I ) » M .NUMY • ETATAB » ZETATAS ) 

M«2 vC 

10 CONTINUE^^^^ 

INITIAL F(XI1)S 

UES2HE1=UEDSTAB(1 ) ^ 

RERS 1 sRERSTAB ( 1 ) ^ 

Rl =RTAB( 1 ) 

DUFDXl =DUDXTAB ( 1 ) 

1F(ZETWTAB(2) >NE«0) DZWOX1 gpZDXTAOn ) 

CAPU 1 *OUEDX I /UES2HE I 
CAPPl «-CAPUl 

IF(ZETWTAB(2) «NE«0) CAPZ1=DZWDX1/( 1 .-ZETAWl ) 
XXL1=XL<1 ) 

COMPUTE INITIAL THETAS 


IFrZETWTABCai -NE.O. JGO TO 11 O 

00 13 I»1*NUMETA C ^ 

13 THETAl ( I >*ZETA ( I ) J 

GO TO 12 

11 DO 201 1=1 .NUMETA 

THETAl ( I )*(ZETA( I l^ZETAWl )/( 1 .-ZETAWl ) 
201 CONTINUE 

12 CONTINUE 




Removed 

M=2 

if(gtab(2).eq.o,)go to it 

Instead of IF(Gl{l-l) ,GT. ,9)M=1 


Removed 

GO TO 16 



17 Gl(l) = 0.0 

16 IF(ZETATAB(2).EQ.O.) GO TO IS 
Removed 


GO TO 19 

18 IF ( ABS ( RHOEROl { I-l ) -1 . ) . GE , . 1 )M=1 

CALL FTLUP ( ETA ( I ) , RHOEROl ( I ) ,M , KUMY , ETATAB , RHOTAB ) 

19 CONTINUE 

ZETAI*n.=ZETWTAB(l) has changed position in program 


Instead of 

IF(ITHETA.EQ. 2 )G0 TO 11 
THETA1(1)=0.0 
DO 200 I=2,I^ETA 

THETAl ( I )=( RHOEROl (l)*(l .-UES2HE1**2-WEDS2HE**2)-ZETAW1+UES2HE1 
1**2*F1 (I )**2+WEDS2HE**2*Gl (I )**2 )/( 1 . -ZETAWl ) 

200 CONTIiJUE 



CO 


c 

c 

c 

c 






c 

c 

c 


c 

r, 

c 


© 


IF 2ETAS ARF GIVE N» INITIAL RHOEROl S MUST 
BE CALCULATED 


ZW«ZETAW1 

IF(ZETWTAB(2 ) .EQ.O. )?W = 0. ^ 

DO 29 I»1,NUMETA 

RHOEROl ( I )= C (1 •-^)»THETA 1 ( I )+ZW-UES2HE1**2*F1 ( I 1**2 4 . 

1-WEDS2HE**2*G! t I 1**2 )/( 1 • “UES2HE 1 ♦*2-WEDS2HE**2 ) 

29 CONTINUE ^ 

DO 21 I*1«NUMETA^ 

RHORHOEI n*l • /RHOEROl < I ) 

21 RH0ER02 ( I )=RH0ER01 ( I ) 

COMPUTE INITIAL VS 


Removed IF(ZETATAB(2) .BQ.O.) GO TO 28 
Instead of ZETAWl 

Removed 28 COR’TIIRJE 


XIBARl = ( 2 *»XI 0 )^HK 2 .»XNBAR ) 

VA M )=(2*XI1 )»»XNBAR /R 1 *»J*RHORHOE ( 1 ) *RERS 1 *VWTAB ( 1 ) /RMUTAB ( 1 ) 

DO 20 I=2»NUMETA^^ 

VA( I )*VA{ I-l )-DELETA ( I-l )*XIBARl (FI ( I )+Fl ( I-l ) ) *XNBAR /(2.*XI0 > 
20 CONTINUE 
YDL( 1 )*0.0 
FNCRERl »RH0ER02 ( I ) 

DO 42 I»2*NUMETA 
FNCRFR2=tRH0ER02 ( I 1 

SUMRERI I )«SUMRER( I-l )+ (FNCRER2+FNCRER1 )/2.*DELETA< I-l ) 

FNCRERl =FNCRER2 

YDL ( I )* (2.*XI 1 )**XNrAR *1 ./RERSI^SUMRERC I ) / ( CAPRS*Rl ♦♦ J*UES2HE 1 ) 
4? CONTINUE 


Removed IF(GEE2.NE.0. )VA(l)=VA(l)*VffiDS2HE/UES2HEl 


calculate INITIAL MS 


VWA=VWTAB ( 1 ) ^ 

DO 44 1*1 *NUMETA 

HOCAPHEI I ) = ( 1 .-^)*ThETA 1 ( I ) -f 2W -UES2HE 1 »»2*F 1 ( I )**2 

1-WEDS2HE**2*G1 ( I )**2 
44 CONTINUE 

CALL CALCM(ZETAW1 tUESZHEl »F1 , WEDS2HE • ZET AF f 
1 SHEfPR* XM8AR1 ,XMC1RC1 ,XMSTAR1 *XMPRIM1 ,PHIR, 

2HSHE* RERSI *Xll *XNBAR .Rl*J, 

3DELETA*RURDRUS* C APRS * SUMRER « AP ♦ BP » CP * ABTAB * FCFT AB • NFCF AB * 

4NMAXGfGl t XXLl > 

40 CONTINUE 


Hei .ovea IF(GEE2.NE .0.)VWA=VVJA*ViEDS2HE/UES21iEl 
Instead of ZETAWl 


C 



(D 


C BEGIN NEW INTERVAL 

C 

1F(XI 1 .GE*X10+30.#DElXUAND. INIT.EO. 1 )G0 TO 46 
GO TO 48 
46 DELXI*DELX10 
INiTsO 
48 CONTINUE 

IFtXI 1 .LT.XITEST )GO TO 35 
DELXI=DELXI*10. 

X1TEST=XITEST*1 0. 

35 XI2=XI1+0ELXI 

IFdCOUNT.EQ.l )XI2=XII 
XIA= fXll+Xl2 )/2* 

ITERATE*! 

DO 30 I*UNUMETA 
F2 ( M=F1 ( I) 

FA ( I )=F1 ( n 
G2 ( I >*G1 ( I ) 

THFTA2( I laTHETAl ( I ) 

THETAA ( I JsTHETAI ( 1 ) 

XMBAR2 ( I )=XMBAR1 ( I ) 

30 CONTINUE 
52 CONTINUE 
C 

C INTERPOLATF FOR F(XI2)?^ 

C 

CALL FTLUP(Xl2»XXL2tl » NUMX * X I TAB * XL ! 

CALL FTLUPIXIA,RA. 1 , NUMX . X I T AQ t RT AB ) 

CALL FTLUP(XIA.RERSA,1 tNUMX . X I T AB . RERSTAB ) 

CALL FTLUP (Xl A ,RMRRMSAt 1 » NUMX • X I TAB « RMUT AB ) 
CALL FTLUPIXIAtVWA. 1 ,NUMX . X I TAB • VWT AB ) 

CALL FTLUP(X!A<RRUUERA«1 »NUMX f XI T AB « RRQUER ) 

CALL FTLUP(XI2»UES2HF2d »NUMXtXI TAB » UECST AB ) 

CALL FTLUP(Xl2tRERS2, 1 tNUMX,XlTAB ,RERSTAB ) 

CALL FTLUP(XI2.R2f 1 »NUMX*XI TAB.RTAB ) 

CALL FTLUPIX I2<DUEDX2. I , NUMX t X 1 TAB . DUDXT AB ) 
IF(ZETWTAB(2 ) •EQ«0« )G0 TO 203 

CALL FTLUP(XI2.ZETAW2« 1 *NUMX*XITA8,ZETWTAB) ^ 
CALL FTLUP(XI2*DZWDX2» 1 • NUMX . X I TAB , DZDXT AB ) 

CAPZ2=0ZWDX2/n .-ZETAWa) ^ 

203 CONTINUE 

UES2HEAS (UES2HE1+UES2HE2 )/2. 

CAPU2*OUEOX2/UES2HE2 

CAPP?=-CAPU2 


00 


Changed position 
Changed position 



0 0 




XIRAR2** (2.»XI2)** (2**XNRAR) 

XlBARA^s (X!0ARI+XIBAR? J/2. 

XIBCPUA= (XIBARl *CAPU1+XIBAR2*CAPU2)/2- 
XIBCPPA= (XTBARl*CAPPl +XIBAR2*CAPP2)/2. 

IF(ZETWTAB(2 )*NE»0. ) XIBCP2A= ( X I BARI -»CAPZ 1 +X I BAR2-»CAPZ2 ) /2 • 
?FnCOt/NT,EQ.l)GO TO ^ 

iterate with SAME XI2 

37 CONTINUE 
7W=ZETAW2 

IFfZETWTAB(2 )#EO.O. )7W=0. 

\t/nO 36 I«lfNUWETA 

RH0ER02 ( I )* M 1 .-ZW>*THETA2( I )+2W-UES2HE2*»2»F2 < 1 )**2 ^ 

1 -WEDS2HE**2*G2 ( I >**2 )/( 1 •-UES2HE2**2-WEDS2HE**2 ) 

RHOeROA ( n = (RHOEROl ( I )+RHOER02 (1 ) >/2* 

RHORHOE (!) = !• /RHOERO? < I ) 

36 CONTINUE 

FNCRERl -RH0ER02 ( I ) 
no 30 I=2*NUMETA 
FNCRER2=Rh0ER02( I ) 

SUMRERI I )=SUMRER( 1-1 )+(FNCRER2+FNCRERl )/2.*DELETA< 1-1) 

FNCRFRl =FNCRER2 

YOL (I ) = (2«*XI2)**XNrAR ♦ 1 • /RERS 2*SUMReR ( I ) / ( CAPRS*R2*» J*UES2HE2 ) 

38 CONTINUE 

calculate ms 


DO 39 I=1*NUMETA 

HOCAPHE( I >*< 1 .-^)*ThETA2{ 1 )+^-UES2HE2**2*F2 ( I )**2 ^ 

I -WE0S2HE*-»2#G2 ( I ) **B 
39 CONTINUE 

CALL CALCM(ZETAW2»UES2HE2,F2 »WEDS2HE.ZETAE» 
t SHE»PRt XMBAR2*XMC IRC2 <XMSTaR2*XMPRIM2«PHIR« 

2HSHE« RERS2«Xt2,XNbAR »R2*J* 

3DELETA .RURDRUS^ C APRS * SUMPER . AP , BP » CP * ABTAB* FCFTA8 . NFCFA0 * 

4NMAXG * G2 > XXL2 ) 

DO 54 I=UNUMETA 

XMBARA ( ! ) c (XMBARl ( I >+XMBAR2 ( I ) )/2. 

IF (GEE2*NE#0. )XMCIRCA I I ) = IXMCIRCl ( I )+XMC IRC2 (I ) > /2 • 

XMSTARA ( I )a (XMSTARl ( I )+XMSTAR2 C 1 ) )/2« 

X14PRIMA n )=r (XMPRIMI < I )+XMPRIM2( 1 ) )/2# 

54 CONTINUE 

NUMDELE = NUMETA- 1 


Removed IF(GEE2.KE.0. )V\^IA=m*WEDS2HE/UES2HEA 


Removed 37 
Instead of ZETAV12 


Instead of ZETA’//2 



DO 50 i*i*numde:le 

XMBARA2 ( M * < XMBARA f I + 1 > +XMBA RA (I ) ) /2 • 

IF (GEE2«NE«0. )XMCIRA2( I ) = fXMCiRCA ( I + l ) + XMCIRCA ( I M/2. 

XWSTRA2(I )*<XMSTARA( T+1 )+XMSTARA( I ) )/2. 

XWPRM A2 ( I) = ( XMPR I MA ( I + 1 ) +XMPR I M A ( 1) ) /2 • 

50 CONTINUE 

calculate FS 

ICHD*1 

FPREV=F2f2) 

55 CALL ABC0GS<NUMDELE*0ELETA*X IBARA.DELXI * FA , VA * XM0ARA2 • ICHD«FI • 

1 X IBCPUA « X I BCPPA •RHOEROA » G 1 * THETA 1 , X 1 8CPZA » THETAA * UES2HEA , XMPRMA2 « 
2F2 * WFDS2HE « G2 ♦ C APG » SmLG * XX ) 

CALL COMPUTE ( NMAXF. FtEST • NUMDELE tOELETA .EPSLONE.C APG « SMLG* I_C^«F2 ) 

update fas and VAS 

DO 70 laE.NUMETA 

FA ( I )» (FI ( I >+F2 ( I ) )/2« 

70 CONTINUE 

DO 73 I«1*NUMDELE 

FA2 { I >» (FA( I+l )+FA ( I ) 1/2* 

73 CONTINUE 

ttmfs=xnbar/xia 

VAd )a (2*XIA)**XNBAR /RA*»J/RhOEROA ( 1 )-*ReRSA»VWA/RMRRMSA 
DO 74 I*2*NUMETA 

VA< I )sVA(I-l )-DELETA(I-I )+X1BARA/(2.*0ELXI )*(F2( I ) +F2 (1-1 ) -FI Cl J 
1-Fl ( I-l ) )-DELETA( I-l )*X 1 BARA+FA2 ( I-l )*TIMES 

74 CONTINUE 

calculate GS 

IF(GEE2.EO.O» )G0 TO 75 ^ 

ICH0*2 ^ 

GPREV=G2(2 ) 

CALL ABCDGS(NUMDELE»0ELETA*X1BARA,DELXI ,FA*VA*XMCIRA2« ichd*fi ♦ 

IX IBCPUA *x I BCPPA *RH0ER0A,GI * THETA l ,X 1BCP2 A , THETAA ♦UES2HEA *XMPRMA2 * 
2F2 *WEDS2HE*G2 * CAPG. SiwLG. XX ) 

CALL COMPUTE ( NMAXG * GTEST * NUMpELE « DELETA * EPSLONE * CAPG * SMLG * 1 CVtP * G2 ) 

calculate thetas 


75 CONTINUE 


CJl 


Instead of rF(GTAB(2).EQ.O. )G0 TO 75 



►1^ 

Oi 



c 

c 

c. 

r 


c 

r 

C 

C 

c 


(D 


ICMO*3 

THFTAPR*THETA2 (3 ) 

CALL ABCDGS (NUMDELE*DELETA»X 1BARA,DELXI tFA • VA * XMSTRA2 » ICHD«Fl ♦ 
1XIBCPUA*XIBCPPA»RH0ER0A*GI ♦THETA! , X IBCPZA ♦ THETAA ♦UES2HEA ♦ XMPRMA2 ♦ 
2F2 ♦ WEDS2HE • G2 ♦ C APG * SmLG • XK ) 

IF (XI I •EQ.XIO.AND* 1 TFRATE.EQ. 1 )NMAXTH=1 . 1»NMAXF 

CALL COMPUTE (NWAXTHfTHTEST.NUMDELE. DELETA » EPSL ONE ♦ CAPG ♦ SMLG . 

I ICHO fTHETA2! 

UPDATE THETAS 

DO 76 1 «1 ,NUMETA 

76 THETAA( I )»(THETA] ( I )+THETA2 ( I ) )/2. 

ITERATE* ITERATE+1 
1F(ITERATE.EQ*2)G0 TO 37 

CONVERGENCE CRITERIA 


IF(ABS( (F2(2 )-FPREV)/FPREV) .LE.EPSLONWIGO TO 90 
GO TO 37 

90 IF(GEE2.EQ*0. )GO TO 1 00 

IF (ABS( (G2(2 )-GPREV)/GPREV) •LE.EPSLONW)GO TO 100 
GO TO 37 

100 IF(A8S{ (THETA2(2)-THFTAPR)/THETAPR) .LE.EPSLONW >GO TO 110 
GO TO 37 
110 CONTINUE 

OUTPUT 


COMPUTE RHOERHO S 
00 80 I=1»NUMETA 

RH0ER02( I )*( ( 1 ♦-2W1*THETA2( I H>ZW-UES2HE2»»2»F2 ( 1 )**Z- 4 

1 WEDS2HE**2*G2 ( I ) **2 ) / U •-UES2HE2**2-WEnS2HF»»2 ) 

RHORHOEf n*I •/RH0ER02( r ) ✓ 

XMMEF ( I ) =F2 ( I ) #SQRT ( RHORHOE t I ) ) 

80 CONTINUE 

FNCRERl SRH0ERO2 ( 1 > 

DO 102 I=2*NUMETA 
FNCRER2»RH0ER02 ( I 1 

SUMREPf I J*SUMRER( 1-1 >+(FNCRER2+FNCRERl )/2.»DELETA ( I-l ) 

FNCRERl =FNCRER2 

YDL ( I )= (2*»X 12 )»*XNbAR *1 . /RERS2*SUMRER ( I ) / ( CAPRS*R2**J»UES2HE2 ) 


Instead of ZETAW2 


Removed IFCGEEa.WE.O. }XMMEG(l)=G2(l )«SQRT(RH0RH0E(T) ) 



102 CONTINUE Z 

IF(XLPR(2).EQ*0.)GO TO 116^ 

IF (XXL2 »GE*XLPR(NPR!NT) )G0 TO 103 
1F(IC0UNT,EQ*1 )G0 TO 103 
GO TO 115 

116 IF(NPRINT.EQ.NSTEPS)GO TO 101 
NPRINTsNPRI NT+ 1 
IFdCOUNT.EQ.l )G0 TO 103 

GO TO 115 

lOi mprint=i 

103 CONTINUE 

DO 91 I »1 fNUMETA 

2ETA( I )=THETA2 ( n*M ^ 

91 CONTINUE 

WR I TE ( 6 f 1 1 1 ) X I 2 » XXL2 

1 1 1 FORMAT! 1H0«3HXI**E20, 8 *4X*4HX/L= ♦E20.8 ) 

MAX=NMAXF+1 0 

GO TO (117*118,119) I VEG 

1 17 WRITEC6* 1 12 ) 

1 12 FORMAT ( 1H0.7X*3HETA • 1 1 X * 3HY/L , 1 2X ♦ 1 HF* 1 3X* 1 HV* 1 1 X ♦ 5HTHET A ♦ 9X , 

1 4HZETA , 0X,8HRHO/RHOE *SX, 5HM/MEF ) 

DO 113 1*1 ,MAX 

WRITE16* 1 14 ) I *ETA ( I ) *YDL ( I ) * F2 ( 1) * V A ( 1 ) ,THETA2 ( I ) ,2ETA ( I ) , 
1RH0RH0E( I ) ,XMMEF( I ) 

1 13 CONTINUE 
GO TO 115 

114 FORMAT! 14,9E14. 4 ) 

118 WRITF!6* 141 ) 

141 format ! 1H0*7X*3HETA * 1 I X , 3HY/L ♦ 1 2X , 1 HF , 1 1 X , 6HEPSDMU • 8X ,5HTHET A * 

I 9X * 4 HZET A , 8X , SHRHO/RhOE * 0X * 5HM/MEF ) 

DO 142 r*l,MAX 

WRITE! 6* 1 14 ) I .ETA ( I ) ,YDL ! I ) .F2 ( I ) ,EPSD«U( I ) .THETA2! 1 ) .ZETA! I ) , 
IRHORHOE! I ).XMMEF( I ) 

142 CONTINUE 
GO TO 1 15 

1 19 WRITE!6. 143) 

143 FORMAT! 1 HO* 7X.3HETA. 1 1X.3HY/L. 1 2X , 1 HF ♦ 1 3X . 1 HG . 1 1 X . 5HTHETA ♦ 9X . 

1 4HZETA.8X. SHRHO/RHOE .8X.5HXMBAR. 1 1X.5HXMMEF) 

DO 144 I si, MAX 

WRITE! 6*114) I ,ETA ! I ) ,YDL (1 ) *F2 ( I ) , G2 (I ) , THETA2 ( 1 ) .ZETA( I ) , 

1 RHORHOE ! I ) . XMBAR2 (I ) » XMMEF ( I ) 

144 CONTINUE 

1 15 CONTINUE 


Removed CALL FTLUP(XI2,XXL,1,NUMX,XITAB,XL) 


Instead of ZETAW2 



update values 


XIBAR2* (2*»XI2)** (2.»XNBAR> 

DO 120 I=1,NUMETA 
FI ( 1 )»F2( I > 

G1 ( T >=G2( I ) 

theta 1 ( I )=THETA2( I ) 

RHOERO 1(1) *RH0ER02 ( I ) 

120 CONTINUE 

UES2HE1 »UES2HF2 
ZETAWl =ZETAW2 
DO 106 I=1*NUMETA 

HDCAPHE ( I ) = ( 1 #->ZW)*THETA1 ( I ) -fZW-UES2HE 1 »»2»F 1 ( I )**2 ^ 

l-WFnS2HE*»2*Gl ( I )**2 
106 CONTINUE 

CALL CALCMIZETAW2* UES2HE1 *F 1 * WEDS2HE » ZETAE* 

I SHE*PR« XMBARl *XMCIRC1 .XMSTARl .XMPRIMl ,PHIR. 

2HSHE* REPS2*Xr2fXN0AR *R2*J. 

3DELETA,RURDRUS* C APRS . SUMRER • AP , BP * CP ♦ ABT AB ♦ FCFT AB ♦ NFCF A0 , 

4NMAXG»G1 *XXL2 ) 

IF ( (COUNT. EQ. 1)G0 TO 104 
IF(XLPR (2 ) .EQ.O. )G0 TO 1 4S 
1F( XXL2 .GE.XLPR(NPRInT > )G0 TO 105 
GO TO 123 

lOB NPPl NT=NPRI NT+1 
GO TO 104 

145 IF (NPRINT.NE. 1 )GO TO 123 
104 CONTINUE 

X I NB= ( 2 .»X I 2 > ♦♦XNBAR 

CFF = PHIR( 1 J/RERS 2 * 2 .»R 2 »*J/XINB* (F2 (2 )'F2 < 1 ) )/DELETA ( 1 )*RuRDRUS 
SUMF = 0.0 

FUNCFl s=F2 ( 1 )♦ ( I .-F2 ( I ) ) 

DO 9? I=2.NMAXF 
FUNCF2=F2 ( 1 )* ( 1 .-F2 ( | ) ) 

SUMF=SUMF+(FUNCF2+FUnCF 1 )/2.*(ETA ( I )-ETA ( I-l ) ) 

FUNCFl =FUNCF2 
92 CONTINUE 

RCURJ=RERS2*CAPRS*UE?2HE2*R2**J 
TSDLF=X INB/RCURJ*SUMF 

GO TO 150 ^ 

SUMG=0.0 

FUNCG 1 -G2 ( 1 ) * ( 1 • -G2 ( I ) ) 
no 94 I=2,NMAXG 
FUNCG2=G2 (I ) ♦ ( 1 .-G2 ( I ) ) 




SUMG»SUMG+(FUNCG2+FUNCG1 )/2.*<ETA ( I )-ETA ( 1 -1 ) ) 

FUNCGl SFUNCG2 
94 CONTINUE 

TSDLG«X I NB/RCURJ#SUMC 

150 CONTINUE 
SM0EL=0.0 
SMTHFsO.O 
FNCDELl*! .“EBn ) 

FNCTHE 1 »F2 < 1 ) ♦ ( 1 • -F2 ( 1 ) ) 

DO 56 I*2*NMAXF 
FNCDEL2*! .-FEI I ) 
c-NCTHE 2®F2 n HK 1 • -F2 ( I M 

SMDEL = SMDEL+ (FNCDEL2+FNC0EL1 )/2.» (YDL ( I )-YDL ( I~1 ) ) 
SMTHE=SMTHE+(FNCTHE2+FNCTHE1 )/2**‘(YDLn )-YOL( I-l ) ) 
FNCDFLl =FNCDEL2 
FNCTHE 1 =FNCTHE2 

56 CONTINUE ^ 

GO TO 151 ^ 

SMDEL=0.0 

SMTHE^O.O 
FNCOELl*! •-G2M ) 

FNCTHE 1 »G2 ( 1 ) # (1 . -G2 M ) ) 

00 57 r»2.NMAXG 
FNC0EL2*! .-G2 ( I ) 

FNCTHE2SG2 ( 1 ) ♦ ( 1 • -G2 M ) ) 

SMDEL»=SMDEL+ (FNCDEL2+FNCDELn/2** (YOL (l)-YDL (I-l)) 
SMTHE=SMTHE+(FNCTHE2+FNCTHE1 )/2**(YDL (I)-YDL M-1)) 
FNCDFLl =FNC0EL2 
FNCTHE 1»FNCTHE2 

57 CONTINUE 
H1SG*SMDEL/SMTHE 

151 CONTINUE 
SUMFRPsO.O 

FNCFRR 1 » 1 . -F2 ( 1 ) *RHORHOE I I ) 

DO 9:^ I=2*NMAXF 

FNCFRR2= 1 • -F2 ( 1 ) *RHORHOE { 1 ) 

SUMFRR = SUMFRR+ CFNCFRR2+FNCFRRI )/2.» ( YDL ( I )~YDL (I-l ) ) 
FNCFRR 1*FNCFRR2 
9-} CONTINUE 

OSDLFsSUMFRR 

GO TO 96 < 

SUMGRR=0*0 

FNCGRR 1 » 1 • -G2 ( 1 ) *RHORHOE ( 1 ) 

DO 97 I»2*NMAXG 


CO 


iienoveQ iiXSF=SMJDL/.”*lTIlE 
Instead of IFCGEES.EQ.O. )G0 TO 151 


•' Instead of rF(GEE2,EQ.O. )G0 TO 96 


FNCGRR2= 1 • -G2 ( I ) ♦RHORHOE ( I ) 

SUMGRRsSUMGRR+ (FNCGRR2+FNCGRRn/2,» ( YOL (l)-YDL (I-l)) 

FNCGRRl =FNCGRR2 
97 CONTINUE 

osdlg=sumgrr 

96 CONTINUE 

REX=CAPRS-»XXL2»UES2Hf2*RERS2**2/RURDRUS 

RDXL=REX/XXL2 

RFTSF=RDXL*T50LF^__ 

REDSF = RDXL*DSDLF^ 

HFDCPHE*! .-UES2HE2**?-WEDS2HE**2 

0BAR=PH1R( 1 )*RURDRUS*R2**J*UES2HE2*CAPRS/(PR*XIN0)*(ZETA{2 ' 

1 -ZETA (1 ) )/DELFTA ( 1) 

ST=O0AR/(RERS2»CAPRS*UES2HE2*(FRO# ( I *-HEDCPHE )+HEOCPHE-ZETA ( i ) ) ) 
WRTTP(6»121 ) 

121 FORMAT! IH0«7X*3HCFF* i OX . 9HTHETA#/LF ♦ 7X ^ 9HDELT A*/LF ♦ 9X ♦ 3HREX ♦ 1 IX t 
19HRETHETA*F,7Xf 9HREDFLTA*F.9X«4HQ0AR. 1 3X * 2HST ) 

WRITE! 6» 1 22 )CFF*TSDLr*DSDLF»REX«RETSFiREDSF fOBAR »ST 

122 FORMAT (8E1 6.4 ) 

WRITF(6* 124 1 

124 FORMAT (1 HO < 7X . 4HH 1 *F i 1 1 X * 3HCCW ) 

WRITE (6* 122 ) HI StCCW 

GO TO 98 ^ 

WRITE(6« 131 ) 

I 31 FORMAT ( IH0*7X.3HCFG< 1 OX ♦ 9HT AUF/T AUG • 6X . 9HTHET A*/LG . 7X , 
19H0ELTA-»/LG»7X«9HRETHETA»G»7X*9HREDELTA*G.9X.4HHI*G) 

WRITE(6« 1 22 ICFG.TAFDTAG. TSDLG.DSOLG.RETSG .REDSG.HI SG 
90 CONTINUE 

WPITE(6.132) 

132 FORMAT nH0*5X.5HYSIDL • 1 1X.5HYS2DL ) 

WRITE!6* 122)YS1DL»YS2DL 
ICOUNTss lCOUNT+1 

123 CONTINUE 

IF<XISTOP.LE.X12lGO TO 130 

XI 1=XI2 

XI0AR1 =XI8AR2 

CAPUl *CAPU2 

CAPPl =CAPP2 

IF(ZETWTAB(2) .NE.O. ) CAPZlsCAPZ2 
, XXL1=XXL2 
GO TO 40 
130 CONTINUE 
STOP 
FNO 


honovpu ...j . 0 . )KhTGG=HDXI * ’Gi)LC*VriLL)S2IiL/UL:r,PhJ' I 


Instead of TF(GEE2.EQ.0. )G0 TO Q& 



SUBROUTINE ABCDGS ( NUMDELE t DELETA • X IBARA * DELX I * FA t VA • CAPMA3 t ICHO« 

I FI «XI0CPUAtXl8CPPA,RHOEROA*Gl tTH.'TTAI *X 1 3CPZA . THETAA .UES2HEA , 

2XMPRM A2 *F2 • WEDS2HE ♦ G? ♦ CAPG ♦ SMLG * XK ) 

D I MENS I ON DELETA { 350 ) « F 1 ( 350 ) ♦ F2 ( 350 ) ♦ FA ( 350 ) « G I ( 350 ) • G2 ( 350 ) * 

1 THETA 1 ( 350 ) * THETA2 ( 350 ) * THETAA ( 350 ) ♦ VA ( 350 ) t CAPMA2 ( 350 ) ♦ CAPA ( 350 ) » 
2CAPB (350) tCAPC(350) , CAPD ( 350 ) ♦ RHOEROA (350 ) .XMPRMA2 (350 ) .CAPG (350) . 
3SMLG ( 350 ) . ZETWTAB (75 ) 

COMMON/ZETW/ZETWTAB 
FCTl =t2.*XK/( ! .+XK ) 

FrT2=2./n . + XK) 
no 10 I=2.NUMDELF 

DELDFL?=2.»DELETA ( I )*DELETA ( 1-1 ) 
yiDFLFA=XTBARA/DFLXI»FA( 1 ) 

VADLDL2=VA ( 1 ) / ( 2 • ♦ ( DFLET A ( I )+DELETA ( I-l ) ) ) 

TMDLDL2= (CAPMA2 ( I ) *Fr T2+CAPMA2 (1-1 )*FCT1 ) /DELOEL2 
CAPA ( I )=VADL0L2 -CAPMa 2C 1 1 *FC T2/DELDEL2 
CAPB ( I ) =X lOELFA+TMOLOLZ 

CAPC ( I )=- (VADLDL2+CAPMA2 ( I-l ) *FC T 1 /DELDEL2 ) 

IF ( ICHO.FQ. 1 )GO TO 20 
IF ( ICHD.FQ.2 )GO TO 3n 
IF( irHn.EQ,3)GO TO 40 

20 CAPD( I )=-CAPA( I )#F1 (1+1 ) + ( X 1 DELF A-TMDLDL2 ) *F 1 ( I l-CAPC ( I )*F1 ( I-l ) 
1-XIRCPUA*FA ( I )#*2 -X IpCPPA*RH0FR0A ( 1 ) 

GO TO 10 

30 CAPO( I )=-CAPA( I )*G1 ( I+l )+ (X 10ELFA-TM0LDL2 )*G1 ( I )-CAPC( I )»G1 (I-l ) 

GO TO 1 O 

40 CAPO ( I ) =;-C APA ( | )*THETA1 ( 1+1 ) + ( X 1 PELF A- TMOLQL2 ) »THETAl ( 1 )~CAPC ( I )* 

1 THFTA 1 ( 1 - i'^UFF2h=-A«-»2/DELDFL2* 

2 (XMPPMA2 ( I )#FCT2* (PI (I+l )^^♦2 + F2 ( 1 + I ) **2-F 1 ( I )**2-F2 ( I )**P) 

3- XMPPMA2( 1-1 )*FCTl*(rl ( I )* + 2 + F2 ( I ) + + 2-F 1 (I-l )**2 -f 2( I-l )*42) ) 

4- \i'FD'' 2HE»*2/OFLOEL2* {XMPRMA2 ( I )+FCT2* (G1 ( l + I )**2 + G2( I + l )*#2 

5- Gl ( I ) *»?-G2 { I ) ) -XMPRMA2 ( I -1 ) #FCT 1 » (G 1 ( I ) ♦♦2 + G2 ( 1 )**2-Gl ( I -1 ) 

/S**2-G2( I-l ) »*2 ) ) 

IF (ZPTWTAR (2 ) , NE . 0» > r APO ( I 1 =CAPD (I) + X!PrP7A»rA(l)»( THFTA A ( 1 )-l , ) 

1'^ rONTINUF 

IF (ZFTWTAR (2 ) ,Nr.O, , and. lOHO.FQ,? )GO TO 41 
IF (7FTWTAB(2 ) .FQ.O. . ANO, IfHO.NF . 1 )GO TO 41 
TA'^’G (2 1 ^-tada ) /rAPR (2 ) 

'.Mt G ( 2 ) =CAPn (2 ) /CAPP (2 ) 

" GO TO 42 

4 I TAPG (2 ) r-rAPA ( 2 )/ (TAPH (2 )+TAPr ( 2 J ) 

SMI 0(2) =CAP0 (2 1 / (C APn ( 2 ) +r Apr (2 ) ) 

4^ rOMTINUr 

on so I -4*NnMDFi r 


Rrr.-rAPM( i ) »r APr ( i ) »rAPG ( i i ) 
r APr, { I ) - -r APA ( i ) /Mrr . 

SHI r, { 1 ) ( r APD ( I ) -rAP<" ( I ) s ( I - 1 ) ) /F)rr, 


cn 


+a JCP ,A*:'A(J ■;T/.A(7 )-3 . ) 



U1 

in:) 


SUBROUT 1 NE CALCR < ZET AW* UEOS2HE *F2 • WEDS2HE • ZETAE * 

1 SHEtpR* XM8AR*XMCrRC*XMSTAR*XW>RlM*PHrR«HSHE* 

2 RERS*XI *XNBAPtR* J, OEUETA,, 

3RURDRUS, CAPRS*SUmRER*AP.BP*CP* ABTAB » FCFTAB *NFCFAB *N«AXG • 
AGE » XXL ) 

Dt MENS I ON F2 < 350 ) ,G2 { 350 ) • XMBAR ( 350 ) *XMC 1 RC ( 350 ) • 

1 XMST AR ( 350 ) t XMPR 1 M ( 350 ) ♦ PH 1 R ( 350 ) * DELET A ( 350 ) « SUMRER < 350 > * 

2ABTAB C20> *FCFTAB (20) ,EMUSDMU (350 ) * ZETWTAB(75) 
COMMON/EPSOMU/EPSOMU (350 ) 

COMMON/TABLE I /RHOERO? ( 350 ) 

COMMON/TABLE2/YDL (350 ) *RHORHOE(350 ) 
common/three/numeta t nmaxf 

COMMON/FEB 1 2 /VW A 

COMMON/FEB 1 1 /PRTT AB ( 20 ) • YDDPRT ( 20 ) ♦ NYP 

COMMON/HDCAPHE/HOCAPhE ( 350 ) • GEE2 

COMMON/MUUSE/ I USEEMU , MPWEMU 

C OMMON / DEC 1 8 / Y CD L 

COMMON/AGA I n/tdl 

C0MM0N/DEC22/YS1DL*Ys2DL 

C OMM ON / J A N 7/ X O 

COMMON/MAY I 1 /A 1 « A2 

COMMON/ JUNE I /SM * SNT 

COMMON/ JULY2/YDEL0 

COMMON/HI S/HIS 

COMMON/ IFTC/IFTC 

COMMON/CCW/CCW 

common/zetw/zetwtab 

COMMON/FBAR/FBARTAB(20)* IFBLU* YDDFB ( 20 ) *NFBY 

common/ I WLDMP/ ! WLOMP 

HRDCPHE = ZETAE~UEDS2HF**2'-WEDS2HE**2 

RURDRUS-SORT ( HRDCPHE/HSHE ) * ( HSHE+SHE ) / ( HRDCPHE+SHE ) 
1*RFRS»HRDCPHE 
2/HSHE 

RUSDRUR*! •/RURDRUS 
SMRERsO.O 

FNCRERI *RH0ER02 ( I ) ^ 

DO 25 1=2*NMAXF ^ ^ 

FNCRER2»RH0ER02 ( I ) 

SMRER*SMRER+(FNCRER2+FNCRER1 )/2.*DELETA< I-I ) 

IF(F2( n.GE..995)G0 TO 23 
FNCRER1*FNCRER2 
25 CONTINUE 
GO TO 23 
24 CONTINUE 



Removed TF{GEE2.NE.0. )G0 TO 2h 




DO 22 I«2*NMAXG 
FNCReR2«RH0ER02 ( I ) 

SMRER=SMRER+(FNCRER2+FNCRER1 )/2*»DELETA( I-l ) 

IFCG2( n *GE.#995)G0 TO 23 
FNCRERI*FNCRER2 

22 CONTINUE 

23 CONTINUE 
NOEL* I 

YDEL«YDL(NOEL) 

DO 60 Ia$*NMAXG 

IF (G2(I ) •GE.*Ol+*99»CCW)GO TO 61 

60 CONTINUE 

61 YS1DL*YDH1) 

DO 70 JJ*I *NMAXG 

IF (G2{ JJ).GE«#5)G0 TO 71 

70 CONTINUE 

71 YDLGsYOL(JJ) 

DO 6? K»JJ*NMAXG 

IF(G2(K).GE« •99+.0l#CC><)GO TO 63 

62 CONTINUE 

63 YS20L*YDL(K) 

IF (YS1DL.LT*0« )YS1DL»0. 

SMDEL«0.0 

SMTHE»0*0 ^ 

FNCDEL1»1 .-F2 ( I > ^ 

FNCTHE1*F2( I )♦( I «-F2 ( 1 )) 

IPI * M-1 

DO 27 I «IP1 «NWAXF 
FNC0EL2*! .-F2C 1 ) 

FNCTHE2»F2( I )■»(! *-F2 (Ml 

SMDEL*SMDEL+ (FNC0EL2+FNC0EL1 )/2*»(YDL<l )-YOL< I-l )> 
SMTHE*SMTHE+(FNCTHE2+FNCTHE1 )/2.*(YDL(1 >-YDL(I-1 )) 
FNC0EL1*FNC0EL2 
FNCTHEl *FNCTHE2 
27 CONTINUE 
60 TO 28 
26 CONTINUE 

FNC0EL1*1*-G2tl ) 

FNCTHEl*G2<i)*(U“G2{l )) 

DO 11 I*2.N«AX6 
FNCDEL2*! .-G2 C I ) 

FNCTHE2*G2< I )•»{! .-GE { I ) ) 


SMDEL»SMDEL+(FNC 0 EL 2 +FNC 0 EL 1 )/ 2 .*(Y 0 L ( I )-YOL (I-l)> 

SMTHE*SMTHE+(FNCTHE 2 +FNCTHEn/ 2 .»(YDL (I)-YDL (I-l)) 


C71 

CO 


Removed TF(GEE2.NE.0, )G0 TO 26 



oi 



FNCOELl*FNCDEL2 

FNCTHE1-FNCTHE2 

11 CONTINUE 
28 CONTINUE 

HtSaSMDEL/SWTHE 

^*FBARMAXaAP+8P*HIS+CP*HIS*HIS 
IF <XXL»NE«XO )GO TO 7? 

YDOO= { 2.*YCDL+TDL)/YDEL0+A3/.4* ( 1 •+ ( -2.*YCDL-. 5*TDL ) /YDELO ) 

72 CONTINUE 

DO 12 I«1.NU«ETA 

PHIR ( I ) =SQRT (HOCAPHE ( I )/ HRDC PHE)* (HRDCPHE+ SHE ) /IHOCAPHE ( I l+SHE) 
EMUSDMUI I )aRUSORUR*RHORHOE ( I )»RERS/PHIR( 1 ) 

12 CONTINUE 

DO 10 I=1.NUMETA 
YOD>=SUMRER ( 1 )/SMRER 

CALL DI SCOT (YODtYDD »YDDPRTfPRTTAB «PRT TAB.- 1 I .NYP.O *PRT ) 


IF( UEQ.l IGO TO 32 

1F( IFBLU.EQ.l )GO TO l3 

YOLDDEL=YDL ( I l/YDEL 

YS2MYS1 * (YS2DL-YS1DL)*SNT 

YOLMYS 1 * YDEL-YS 1 DL 

IF( IFTC.GE. 1 )GO TO 52 

IF(CCW •GE.«85)G0 TO 52 

IF<YS2MYS1-»A2/A3.GE« YDLMYSl )G0 TO 52 

A1 YCDLsAl*YCOL 

A2YS2M1 sA2*YS2WYSl 

TW0=(.8-A1 )*YCDL/*4*YDEL/YDEL0 

AMADTMY= < A 1 YCDL-*A2YS2M 1 )/(TWO-YDLG ) 

YDOOYDL=YDDO*YDEL 

IF ( YS2MYS1 .LE.YCDL*A1 /A2 )GO TO 50 


IF (YS2MYS1 .GT.YCDL4A1/A2 IGO TO 51 
50 YCOLDOL=YCDL/YDEL 


© < 


0NE=Al*YC0L/.4 

IF(YDL( 1 WLE«0NE)FBAR**4*YDLD0EL 

IF (ONE.LE.YDL ( I ) •ANO.YDL < I ) • LE .TUO ) FBAR- A1»YCDLD0L 
IF(TW0.LE#YDL< I l.AND.YOL ( I ) .LE. YDLG 1 FBARsYDLDDEL*AMADTMY+A 1 YCDL/ 

1 ydel-two*amadtmy/yoel 

IF(YDLG.LE.YDL( I >.AND»YDL( I ) .LE* YDDOYDL )FBAR=A2YS2M1 /YDEL+ (YDL ( I )- 
1 YDLG )/YDEL* ( A3* < YDEL-YS 1 DL ) --AE^YSEM YS 1 ) / ( YDOOYDL-YDLG ) 

IFfYDL( I ) .GT.YDD0YDL)FBARaA3*( 1 •-YS1DL/YDEL ) 

GO TO 53 

51 0NE=A2*YS2MYS1/«4 

IF(YDL( I )«LT.ONE)FBAR=«4*YDLDOEL 

|F(0NE*LE.YDL( I )*AND«YDLn ) . LE . YDLG ) FBAR* A2*YS2MYS I /YDEL 



IF(YDLG«LE«YDL( I ) .ANn«YDL( I ) .LE. YDDOYDL )FBAR=A2YS2M1 /YOEL+ (YDL( I )- 
1 YDLG ) /YDEL* t A3* ( YDEL-YS 1 OL )“A2*YS2M YSl ) / ( YDDOYDL-YOLG ) 

IFCYDLC I ) .GT.YDDOYDL)FBAR= A3* (I •-YS1DL/YDEL) 

GO TO 53 
52 IFTC=IFTC+1 

2007 IF(YDL{ I ) .LE* . 1 ♦YDEL)FBAR* .4*Y0LDDEL 

IF(YDLD0EL*GT« *1 .AND«YDLD0EL*LE# .3 )FBAR* • 04+ (YDLDOEL-. 1 J/.2 
l*(FBARMAX-.04 ) 

IF(YDLDDEL«GT, ,3 JFBARaFBARMAX 
i 53 CONTINUE- 
^ GO TO 14 

13 CALL DISCOT IYDD*YDD« Y0DFB*FBARTAB.FBARTAB*-1 1 *NFBY.O*FBAR) 

14 CONTINUE 

TF( I .EO*NUMETA)GO TO 32^ 

EPSDMUI I )=(2**XI )**XNBAR/R**J* FBAR**2*EMUSOMU ( I )*RH0RH0E ( I) **2* 
1SMRER*»2 

1*ABS( <F2( I+l )-F2( I-l ) )/(DELETA( I J+OELETA (1-1 ) M 
GO TO 34 
29 CONTINUE 

EPSDMUt I )»(2«4X! )**XNBAR/R**J*FBAR**2*EMUSDMU( I ) *RH0«H0E ( I )»*2 
1*SMRFR**2 

l*SORT(ABS( ( (F2( 1+1 )-F2( I-l ) >/(DELETA ( I )+DELETA( I-l ) ) 1**2 
1+ ( (G2 ( I+l )-G2 ( I-l ) )/ (DELETA ( I l+OELETA (I-l ) ) )**2 
1 * ( WEDS2HE/UEDS2HE ) **? ) ) 

34 CONTINUE 

IF( IWLDMP*EO*l ♦AND. I .NE.2)GO TO 15 

TF(I*NE.2)GO TO 33^ 

FOCFD2» ( ( RHORHOE ( 1 )>VW A ) / ( PH I R ( 1 ) /RERS* ( 2 .*R**J ) / ( (2**XI )**XNBAR)* 

1 F2 ( 2 ) /DELETA ( 1 ) *RURDRUS ) ) *2 • 

GO TO 41 

40 FDCFD2=( (RHORHOE! 1 ) *VWA ) / ( PH I R ( 1 ) /RERS* ( 2 • *R**J ) / ( (2«*XI )**XNBAR)* 
1SQRT( (F2(2)/DELETA( 1 ) )**2+ (G2 (2 )/OELETA ( 1 )) **2* ( WEDS2HE/UEDS2HE ) ** 
22)»RURDRUS) )*2. 

2*UEDS2HE/S0RT ( WEDS2HF»*2+UEDS2HE**2 ) 

41 CONTINUE 

CALL FTLUP(F0CFD2*AB,1 #NFCFAB*FCFTAB.ABTAB) y 

15 CONTINUE 

IFMWLDMP,EQ*1 )G0 TO 16 

ADL= AB*S0RT( (2«*XI ) **XNBAR ) / (CAPRS*SQRT ( fmusdmu < 1 ) )*RH0RH0E ( 1 ) 
1*RERS 

1*UEDS2HE*SQRT(R**J)*sQRT(ABS(F2(2) ) /DELETA (1 ) ) ) 

GO TO 36 

16 A0L=AB»SQRT( (2.*XI J **XNBAR )*SQRT (EMUSDMUI 1 ) )/(CAPRS*EMUSOMU( I ) 

1 *SORT ( RHORHOE ( 1 )) *SORT (RHORHOE ( I ) ) *RERS 


Removed IF(GEE2.NE.0. )G0 TO 29 


Removed IF(GEE2,NE.0, )G0 TO UO 


Removed IF(GEE2,NE,0. )G0 TO 35 



0 0 © 


tn 

a> 


{ 


l*UEOS2HE*SOPT(R**J)*saRT(A0S(F2(2) )/DCLETA(l ) ) ) 

GO TO 36 

35 CONTINUE 

ADL* AB*SQRT« <2«*XI ) **XNBAR ) / <CAPRS*SQRT (EWUSDMU< 1 ) )*RMORHOE ( 1 ) 

1 #RERS 

1 *UEOS2HE*SQRT ( R**J ) 

1 * ( ( F2 ( 2 ) /OELETA ( 1 ) ) *#2+ C G2 ( 2 > /DEUETA ( I M ♦♦2* ( WEOS2HE/XrEDS2HE ) **Z ) 
1*».25) 

36 CONTINUE 


33 EPSDMU ( I ) *EPSO«U < 1 ) ♦ { 1 .-EXP < -YOL 
IF <EPSDMU{ 1 ).LT.O. )EPSOMU( I )*0. 
GO TO 31 

32 EPSDMUC I )s0.0 
31 CONTINUE 

IFdUSEEWU.EQ.O )EPSOWU< I )»0. 
XW8AR( I >*PHIRr I .+EPSDWU( I > ) 
XMCIRCt 1 )« PHIR< 1 )* ( 1 ./SW-fEPSOWU< 
XMST AR { I J »PH I R ( I) /PR* ( 1 • +EPSDMU ( 
ZW»ZETAW 

IF(ZETWTA8(2 ).EQ.O. )ZW»0. 

1F( WPWEMU.EQ.l )GO TO 37 
XMPRIMI I )«PHIR( I )/PR*( 1 .-PR)/( 1 . 
GO TO 10 

37 XMPRIMI I )*PHIR( I ) /PR» (i .-pp)/< i . 
IRT ) 

1 / n .-PR j ) 

10 CONTINUE 
RFTURN 
END 


( I )/A0L) )**2 


l/SNT) ^ 

)*PR/PRT ) 

^ 

■ZW )* ( 1 .+EPSDMU( I ) *PR/PRT* ( 1 *-P 


Instead of XMCIRC(l)=XMBAR(l) 


Instead of ZETAW 
Instead of ZETAW 



u u u 


{ 


subroutine compute (NmAX. test tNUMDELE* DELETA *EPSLONE,CAPG *SML<3 t 
1 1CH0*DUM2) 

D I MENS I ON CAPG ( 350 ) . SMLG ( 350 ) • DUM2 ( 350 ) . DELET A ( 350 ) • ZETWTAB(75) 

COMMON/CCW/CCW 

C0MM0N/ZETW/2ETWTAB 

C0MM0N/ZETAW2/ZETAW2 


test-finding edge 


KMAX=NMAX-10 

DO 10 JJ=KMAX*NUMDELE 

TSTV AL =EPSL0NE»DELET A ( JJ ) 

CKVAL = TEST*U .0-CAPG( JJ) )-SMLG( JJ) 

IVAL=*JJ+4 

IF(ABSCCKVAL).LE.ABS(TSTVAL) )GO TO 20 
10 CONTINUE 
20 NMAX-IVAL 

C 

C COMPUTE DUM2S 

C 

NBACK=-<NMAX-1 ) 

M2*-2 

DO 30 NF=NBACK,M2 
KF«IABS(NF) 

DUM2(KF)=CAPG(KF)*DUm2(<F+1 )+SMLGIKF1 
30 CONTINUE 

t IF ( ICHD.EQ.2 )DUM2 ( 1 >=CCWxDUM2 (2 ) 

IFt2ETWTAB(2).EQ.0*.AND« 1CH0.EQ*3)DUM2(1 > «ZETAW2»OUM2 t 2 ) 

RETURN 

END 
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X, cm 



(a) Velocity profiles, x < 11.7 cm (4.6 in.). 



(b) G profiles and lines of constant G. x < 11.7 cm (4.6 in.). 



(c) Mixing-length distributions, x < 11.7 cm (4.6 in.). 


Figure 2.- Theoretical predictions for velocity and concentration profiles and mixing- 
length distributions for s = 1.73 mm (0.068 in.) rearward inclined step slot configu- 
ration of reference 21. Experimental data for velocity profiles are shown for compari- 
son. aj = 0.14 used for all solutions; y^ j and y^, jj evaluated at G = 0.99 and 

G = 0.01, respectively; t = 0.025 mm (0.001 in.) used in calculations. 
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(a) X < 12.2 cm (4.8 in.). 

Figure 3.- Comparison of predicted and experimental velocity profiles from data 
of reference 21 for rearward inclined step slot with s = 0.038 cm (0.015 in.). 
aj=0.14; am -0.09; Ngc,T = 0.8. 
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(b) X > 12.7 cm (5.0 in.). 


Figure 3.- Concluded. 





Figure 4.- Calculated distribution of skin friction for aj = 0.14, a^ = 0.09, and Ng^ ^ = 0.8 and compari- 
son with experimental data from reference 21. Flagged symbols at x = 50.8 cm (20 in.) are skin-friction 
balance data and unflagged symbols are from measurements of 0 on solid flat plate without a slot. 


F or C 

Figure 5.- Initial profiles at Xq based on data from reference 23. 
Mg = 6; Tt,e ~ 444° K (800° R); Ttj = 294° K (530° R); 
s = 0.25 mm (0.01 in.); X = 0,068. 
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Figure 6.- Comparison of computed heat-transfer distribution with experimental 

data of reference 22. 



Figure 7.- Comparisons of calculated values of effectiveness with experimental 

data from reference 22. 
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X = (Ax/s) 

Figure 9.- Comparison of computed heat-transfer distribution with 
experimental data of reference 24. 
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Figure 11.- Comparison of predicted effectiveness with experimental data for 
Me = 6.0 and X = 0.06 to 1.6 (ref. 29) based on direct measurements of 


recovery temperature for s = 0.15 to 1.12 cm (0.06 to 0.44 in.). Values of 
input parameters in theory were X = 0.065, Me = 6.0, aj = 0.14, 
s = 4.78 mm (0.188 in.), a^^ = 0.09, and Nge,T ~ 



Figure 12.- Comparison of Cf prediction from finite-difference solution with 

experimental values obtained by Cary from skin-friction balance. Flagged 
symbols are measured values without slot, s = 4.78 mm (0.188 in.); 
t = 1.57 mm (0.062 in.); X = 0.065; 6 q = 50.8 mm (2.0 in.); 

Tt,j/T,,e = 0.63. 
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